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results are in reasonable agreement with the limited experimental data available. 
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technique Integrates the horizontally directed momentum equation utilizing empirical 
Information about the flow field behind the backward- facing step. The recirculation 
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The results predict velocity profiles in sufficient detail that the presence of the 
corner eddy in the region of negative surface pressure gradient is evident. The 
magnitude of the reversed flow velocity in the recirculation eddy has been found to 
agree with that found from experiments. Also, a surface eddy viscosity distribution 
has been an outgrowth of the method which realistically follows the magnitude of the 
surface pressure gradient distribution as found experimentally. 
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CHAPTER I 


INTRODUCTION AND STATEMENT OF THE PROBLEM 

The study of atmospheric flow over buildings or 
other man-made surface obstructions has gained momentxim in 
recent years due to the possibility of the operation of 
short tcdce-off and landing (STOL) aircraft in addition to 
helicopters from within the heart of a city. Unsteady flow 
phenomena such as induced vortex fields, cross-winds, and 
separated flows near these obstructions pose a danger to the 
operation of low-speed aircraft, especially in a strong 
turbulent wind. In fact, some accidents involving small 
aircraft while taking off from or landing in around airport 
terminal buildings have been traced to these flow phenomena. 
In structural design of buildings or bridges the civil 
engineer has always been concerned about the dynamical loads 
induced by these complex flow phenomena on such structures, 
for the lack of understanding of which, a conservative 
factor of safety has been used in their design. Even then, 
occasional incidents of glass panes being shattered out of 
buildings in strong winds have been reported. Aero-elastic 
phenomena such as buffeting and stalling arise primarily out 
of the separation phenomenon. Vibration problems involved 
with the collapse of the Tacoma Narrows bridge stemmed from 
these flow phenomena in a strong turbulent wind. There have 



also been reports of severely unpleasant situations around 
downtown shopping centers due to strong recirculating flow 
behind the buildings. 

These flow phenomena are so complex that they have 
to be dealt with separately. In this study the problem of 
flow separation behind an obstruction, its reattachment to 
the ground, and the recirculating flow in the cavity or the 
wake zone thus formed has been investigated and solved by 
considering the simple configuration of a rearward- facing 
step. Some experimental data are available for such an 
obstruction and permit computational solutions to be 
compared with experimental findings. To simplify the 
problem further only two-dimensional flow is considered. 

The problem has been attacked two ways. While the 
first method is the solution of general Navier-Stokes 
equations as applicable to turbulent atmospheric flow, the 
second approach is an integral technique to predict the 
velocity profiles in the recirculation region by making use 
of the experimentally known pressure fields. From the 
former approach vorticity and stream function are obtained 
as a solution without taking recourse to the knowledge of 
pressure. However, once the solution is obtained, pressure 
field can readily be achieved. 
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CHAPTER II 


THEORETICAL REVIEW 

In this chapter major theories that have been put 
forward to predict the separation phenomenon of laminar and 
turbulent flows over an obstacle are discussed. Results 
from experiments conducted by various authors are reviewed 
to furnish appropriate empirical data in terms of constants 
and parameters which are used in both approaches to the 
numerical solution of the problem presented herein. 

This chapter is divided into three sections. In 
the first section an introductory paragraph about the 
phenomenon of separation is given. In the second section 
the mechanics of the flow region between separation and 
reattachment and that downstream of reattachment is reviewed 
In the third section various analytical models as proposed 
by different authors to predict the parameters associated 
with the laminar and turbulent separation, for example, 
velocity, eddy viscosity and turbulence intensity in the 
free shear layer, and base pressure are discussed. 

I. SEPARATION 

The phenomenon of separation occurs due to the 
presence of either viscosity and a pressure gradient or of 
an abrupt change of geometry. In the latter case, on most 


3 



structures, it occurs from the sharp edges. Separated flows 
are formed, for example, upstream of a forward- facing step, 
downstream of a backward- facing step, within a cutout in a 
body surface, and on the upper surface of an airfoil at high 
angle of attack. 

When a two-dimensional flow separates due to 
viscosity in an adverse pressure gradient, retardation of 
the fluid close to the surface causes a rapid thickening of 
the boundary layer. The point where the velocity gradient 
at the surface in the direction normal to the wall decreases 
to zero is defined as the point of separation. Downstream 
of this point flow reversal takes place near the surface 
(Figure 2-1). At the separation point the shear stress is 
reduced to zero. Experimentally, the occurrence of zero 
shear stress and the flow reversal near the surface has been 
found for strictly a steady, two-dimensional, incompressible 
laminar flow. When turbulence and three-dimensional effects, 
as in the case of a building, are introduced, this 
phenomenon is not well defined. 

Rearward- facing step- flow separation is charac- 
terized by a relatively small angle of incidence between the 
streamline and the body at the point of separation and at 
the poj nt of reattachment; but in the case of bluff-body- 
flow separation, the streamline at separation and 
reattachment points is inclined almost perpendicularly to 
the body and the reattaching surface. The step-flow and 
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(a) Flow past a body with separation 
(S = point of separation) 



(b) Shape of streamlines near point 
of separation 



(c) Velocity distribution near the 
point of separation 
(PI - point of inflection) 


Figure 2-1. Separation of the boundary layer [22] . 
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bluff -body- flow types of separation arise due to an abrupt 
geometrical change in configuration of the body. Separation 
bubbles on airfoils occur due to the dynamical interactions 
of viscosity and an adverse pressure gradient which sets up 
recirculating motion on the airfoil surface. 

The flow separation depends heavily on whether the 
flow is laminar, turbulent, or transitional. In subsonic 
flow pure laminar separated flow does not hold any practical 
importance [1];^ but for compressible flow, and especially 
at higher Mach number, laminar separation is important, 
since it is very stable and may persist to high Reynolds 
numbers. In this study, however, only low Mach number 
incompressible flow is considered, 

II. MECHANICS OF THE FLOW REGION 

Flow Between Separation and Reattachment 

Shear and turbulence in the incident flow are two 
important features of flow separation around a surface 
obstruction o The distorted flow over a block geometry 
building consists of a displacement zone, a wake region 
which encloses the rear separation bubble, and an upstream 
separation bubble (Figure 2-2) . Shear generates a swirling 
flow in the wake or cavity zone, which can be interpreted as 
the accumulation of vortex lines in the wake region by the 


^N\imbers in brackets refer to similarly numbered 
references in the bibliography. 
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Upatrcan nopariitlon Rear eoparailon bubble 

bubble or down waen or cavity none 

Bone 



Reattaoheect flow r.oiie 


Figure 2-2. Definition of flow zones near a sharp 
edged building [2], 
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passing flow. Turbulence thickens the shear layer which 
encloses the wake region. For some distance over the 
obstruction, and in the wake zone, turbulence is very large. 
The fluctuating component can be as high as 40% of the 
steady component [2] and is generally largest on the edge of 
the shear layer. In most cases the vortices of the sepa- 
rated flow are unsteady, and their experimental study is 
difficult. 

Wind tunnel studies of separation over rearward- 
facing steps and bluff geometries [1 through 11] provide 
some insight into the physical aspects of flow separation 
and recirculating flow formation in the cavity zone. When a 
two-dimensional turbulent flow separates from 'the sharp 
corner of a back step, a shear layer with high vorticity and 
low static pressure is formed which spreads linearly down- 
stream (for steady, two-dimensional laminar flow the 
spreading is parabolic) Momentum diffuses from this 
turbulent mixing layer into the cavity zone which sets the 
wake fluid into motion, and thus the sharp velocity 
discontinuity at the wake boundary is smoothed out. In the 
reattachment zone the pressure increase arising out of the 
compression creates a steep pressure gradient near the 
surface such that part of the flow near the surface returns 
upstream to feed the recirculation zone. Due to this 
entrainment process the velocity and pressure variations are 
large, and the wake boundary is not well defined. For 
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theoretical analyses, however, this boundary, also called 
the dividing mean streamline, is chosen as the streamline, 
ij; = 0, This line separates the mixing layer, which 
comprises the original and part of the new shear layer, and 
the outer (undisturbed) flow from the recirculation region 
(Figure 2-3) where the fluid recirculates as a large eddy. 
The corner eddy is formed as a result of a shear-layer 
separation which is induced by the reverse flow now 
approaching the step as a forward-facing step. Alternately, 
for theoretical analyses the steady state shape of the 
recirculation bubble may be assumed to arise from the 
balance of the rate of entrainment from the bubble into the 
turbulent mixing layer along the bubble boundary and the 
rate of reversal of fluid back into the bubble [2] . 

Tani, et al ^ [8], used the model shown in Figure 
2-4 for turbulent flow separation behind a back step, 
measuring surface pressure coefficients averaged over a set 
of measurements , 


C 


P 


P-Pr 

; 2 " » 

1/2 pu; 


where p and p^ are the local and undisturbed static 
pressures, respectively. These are plotted versus the 
nondimens ional distance, x/h, in Figure 2-5, where x is the 
horizontal distance downstream from .the step of height h. 
Except for steps of very small heights, the base pressure 
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Figure 2-3. General behavior of typical reattaching 
flow [4] . 
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Figure 2-4. Separated flow model over a backward- 
facing step [8] . 
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Figure 2-5. Pressure distribution on the step face 
and the bottom face [8], 



was found to be insensitive to the step height and to the 
initial boxindary layer thickness, the reason apparently 
being that the cavity flow is chiefly maintained by the 
turbulent shear stress, which is almost independent of the 
step height and the approaching boundary layer. It was seen 
(Figure 2-6) that turbulence and shear stress increase down- 
stream in the mixing region and that the distribution across 
the mixing region of mean velocity, turbulence intensity, 
and shear stress is also insensitive to the step height and 
initial boundary layer thickness. The reattachment length 
was found to be equal to about seven step heights downstream 
of the step. 

Experiments conducted by Bradshaw and Wong [4] 
show that for a wind tunnel speed of 25 meters per second, 
the recirculation zone is about six step heights long for a 
step height of 2.5 cms. It was found that after reattach- 
ment the turbulent shear layer splits up, part of which 
proceeds downstream with the eddy length scale considerably 
reduced, and the other part of which turns back upstream 
toward the recirculating region to supply the entrainment. 
Large eddies in the shear layer are virtually split in two 
if the fraction of the shear- layer mass flow that is 
deflected upstream at reattachment is appreciable, which is 
the case with a large initial boundary layer thickness and 
hence the atmospheric boundary layer. As a consequence, 
the turbulence structure is drastically different from that 


12 



h = 2 cm 

= 28 m/sec 




K 


Figure 2-6. Distribution of streamwise mean 

velocity, turbulence intensity, and turbulent 
shear stress behind a backward- facing step [8]. 
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found in any conventional shear flow. It had been argued 

earlier that large eddies are deflected alternately upstream 

and downstream rather than split as proposed by Bradshaw and 

Wong [4] , This would lead to a stronger unsteadiness in the 

wake than is normally found. At reattachment the thickness 

of the new shear layer is of the order of a step height 

(Figure 2-3, page 10). Curvature of the mean streamlines 

and a rapid distortion due to the pressure field near 

reattachment are the factors causing the initial shear layer 

to split roughly in half at reattachment (Figure 2-3) . 

As measured by Tani , et al , [8], and Mueller and 

Robertson [9], the shear stress in the free shear layer is 

2 

much higher than the value of 0.01 pUj^ (U^ is the free 

stream velocity) found in a plane mixing layer. This is 

true even though the curvature of the streamlines associated 

with a free shear layer tends to decrease the shear stress 

and turbulence intensity. The reason for higher stress i^^ 

that the effective velocity difference across the shear 

layer is more than due to the reversed flow in the 

separated region. Although the reverse flow velocity does 

not seem to exceed 0.2 , the shear stress exceeds 

2 

(1 + 0.2) times the plane shear layer value. Shear stress 
in the reversed flow region is not negligible [7, 10, 11] , 
and it seems that the shear stress in the shear layer is 
increased by the "feedback" or the re-entrainment of 
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stress-bearing fluid from the separated-f low region. Such a 

mechanism needs further experimental investigation. 

9U 

The mean velocity gradient, , on a given stream- 
line will be nearly the same before and after the region of 
rapid distortion at reattachment [4] , that is to say, the 
mean vorticity is nearly conserved along a strecimline; but 
despite this the decrease in Reynolds stress, -u* v* , 
accounts for a sudden drop in the value of maximum shear 
stress near reattachment (Figure 2-7), Evidently large 
changes in turbulence structure occur when the shear layer 
bifurcates at reattachment. 

As has been reported by many authors, for a thick 
initial boundary layer the maximum shear stress occurs near 
the dividing streamline throughout most of the separated 
shear layer except near reattachment where the Reynolds 
stress, m ' V ' , decreases rapidly to a value of zero. The 
turbulence intensity on the dividing streamline also 
decreases, but less rapidly. Although the measurements of 
shear stress near reattachment are not completely reliable 
[4] , it is fairly well known that the surface shear stress 
rises rapidly after reattachment. 

Flow Downstream of Reattachment 

On the basis of experimental evidence, Bradshaw 
and Wong [4] define and classify the perturbation in the 
initially thin shear layer with the value of the parameter. 
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Figure 2-7. Maximum shear stress in shear 
layer [ 4 ] . 
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1i/6q (step-height/boundary layer thickness measured at the 
separation point) , as follows: 

1. Weak perturbation : ^76^ << 1, 

2. Strong perturbation : 1^/6^ = 0(1), 

3. Overwhelming perturbation; 

(This parameter has also been found useful for classifying 
the flow downstream of reattachment.) An example of an 
overwhelming perturbation is the mutation of a boundary 
layer to a wake or a mixing layer. The flow of a thin 
boundary layer over a downstream- facing step involves two 
such perturbations, boundary layer to mixing layer and 
mixing layer back to boundary layer. The first perturbation 
can be ignored, and the flow can be treated as if a fully 
developed mixing layer appeared at the separation point. 

Since the atmospheric shear layer is generally large 
relative to any step height, that is, h/6Q << 1, a building 
in atmosphere could be classified as a weak perturbation. 

An equilibrium boundary layer is defined as one 
for which the turbulence structure is not strongly disturbed. 
The relaxation of the boundary layer after reattachment, a 
process where the disturbed boundary layer returns to its 
equilibrium state, has been found to be very slow [4, 12, 

15], It is generally assumed that the distribution of the 
average velocity can be described by universal laws applied 
to two layers in the turbulent boundary layer, the "Law of 
the Wall" and the "Law of the Wake." If this average 


17 



velocity distribution follows these universal laws, it will 
be assumed that the turbulent boundary layer is in an 
equilibrium state [12], 

The "Law of the Wall" and the "Law of the Wake," 
for instance, as defined by Coles [13] are written as 
follows : 



, zu 

In 

— in 

K V 


+ c + 


2tt 

K 


- • 2 /TT Z» 

sin (jj) 


u : average velocity with respect to z, 

z ; perpendicular distance from the wall, 
u^: /t^/p , shear stress velocity, 

V ; kinematic viscosity, 

TT : a function of pressure, 

K ; 0,4 (von Karman*s constant), 

Townsend [14] has an alternate approach to the 
definition of an equilibrium turbulent boundary layer. 
Townsend defines the state of equilibrium to exist between 
the production of Reynolds shear stress and the dissipation 
of turbulence energy which introduces gradient terms such as 
4^ and -4^ in the analysis. Obviously, the law of the wall 
and wake is a simpler approach. 

The characteristic lengths that describe the 
average velocity in a turbulent boundary layer are v/u^ in 
the range of the Law of the Wall and 6 in the range of the 
Law of the Wake, It is reasonable to assume that the 
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relaxation process could as well be characterized by any two 
lengths [13]. However, in the author's belief [12], 
additional parameters seem necessary to determine the 
process of relaxation completely. Once experimentally 
determined, they could be incorporated into the calculation 
procedures. However, these parameters are yet unknown. 

The relaxation of the shear layer after reattach- 
ment to the ordinary boundary layer has been found to be 
very slow and non-monotonic [4]. Earlier reports suggested 
that the relaxation was monotonic [8, 10]. However, the 
results of [4] show that the surface shear stress does not 
return monotonically to the constant pressure equilibrium 
value and that it takes more than 50 step heights to return 
to its equilibrium value. These results correspond to an 
almost overwhelming perturbation, [^/Sq = 0(10)]. But even 
for just a strong perturbation, [h/S^ = 0(1)], Wauschkuhn 
and Ram [12] found that more than 20 step heights downs trecun 
were required to reach the equilibritim state. The results 
of Wu, et al . [15] , from experiments to investigate different 
features of flow separation over a back step (a strong 
perturbation) also indicate that this relaxation distance is 
very large; and it was observed that static pressure takes a 

shorter distance to return to its equilibrium value than the 
velocity in the shear layer. These results indicate that 
the earlier belief which considered the mean velocity as 



following the logarithmic law of the wall close to the 
surface in the relaxation region is wrong. 

The recovery of the boundary layer after reattach- 
ment depends appreciably on the obstacle shape, even when 
the separation occurs from a sharp edge. In other words, 
obstacle shape affects the strength of the recirculation in 
the cavity zone. Relaxation behavior of the boundary layer 
after experiencing an adverse pressure gradient, for example, 
on a flat plate, which is monotonic in nature, is different 
from that of the boundary layer after reattachment, in the 
case of a step- flow, the primary reason being the rapid 
distortion of the flow in the reattachment region. The 
relaxing boundary layer in the case of a strong perturbation 
experiences more severe disturbance than in that of an 
overwhelming perturbation, although recovery from the former 
is quicker for a given initial boundary layer thickness, 6^. 

After reattachment, the turbulent length scales, 
especially the dissipation length parameter, L (for ordinary 
boundary layer, L = kz) , are almost independent of z except 
for a sudden drop to a value of zero at the surface [4] . 

This reduces the velocity gradient and hence the velocity in 
the inner layer below the value predicted by the logarithmic 
law. 
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which demonstrates that the Law of the Wall is inapplicable 
under these conditions. Alternatively, the local-equilibrixim 
form of the mixing- length formula [14] , 

8z KZ 

gives even higher values of the velocity gradient. The 
failure of this formula can be traced to the fact that 
turbulence is not in local (energy) equilibrium, but changes 
rapidly in the streamwise direction, and that the length 
scale of turbulence is not proportional to z , but increases 
much more rapidly with z near the wall. 

The outer layer, defined as one which retains the 
characteristics of the mixing layer, may take longer to 
return to the normal boundary layer state, since the outer- 
layer eddies, being larger, have longer lifetime than the 
inner-layer eddies. The effects of rapid distortion in the 
reattachment region propagate to the outer layer and change 
its mixing-length characteristics. Due to bifurcation of 
the mixing layer in the reattachment zone, the central 
region of the mixing layer comes into close contact with the 
surface; and hence there is a sudden jump in the value of 
apparent mixing length or true length scale of turbulence. 


L 


(-iPT*') 


3/2 


dissipation rate 


above the local equilibrium value with increasing z. 
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Mueller and Robertson [10] conducted experiments 
on wedge-type roughness elements in air for the element- 
height Reynolds numbers such that the shear layer was 
turbulent downstream of the element. Turbulence was kept 
low by using mesh screens. The reattachment point was 
determined as the location where the shear stress value was 
zero, which was roughly seven element heights downstream. 

The mean velocity profiles behind the elements are shown in 
Figure 2-8, Boundary layer thickness is seen to increase 
slowly with distance downstream. Local skin friction 
coefficients, , were obtained from the profiles near the 
wall from the Law of the Wall given below; 

C, zu , C, 

^ = 5-6 (-^) [1 + log(-^) + 4 log(-^)] 

The results are shown in Figure 2-9. With the assumptions 
of a fully developed mixing region for the largest element, 
shear stress at the edge of the separation bubble, the 
dividing streamline, was calculated from the error function 
curve [16], The calculated and the measured shear-stress 
profiles are compared in Figure 2-10. The width of this 
mixing region was found to increase linearly with distance 
from the roughness element. 

The authors [10] classify the turbulence structure 
behind an obstacle into two stages. The first, one of 
excitation, from separation to the reattachment and the 
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Figure 2-9. Dovmstream redevelopment of 
mean flow [10] , 
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Figure 2-10. Turbulent shear-stress profile 

at x/k = 4.5 behind 0.7 5- inch element [10]. 
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second, one of decay, after reattachment to the equilibriiim 
state downstream. These two regions are clearly defined in 
Figure 2-11. From Figure 2-12 it can be seen that the shear 
stress has an initially high value in the excitation region. 
It decays downstream to the normal boundary layer value. 

For smaller elements there is a different maximum shear 
stress curve for the excitation region which indicates that 
the mixing layer did not reach the fully developed stage. 
Also, interference from the base could affect the shear 
stress curve for smaller elements. For the smallest element 
the boundary layer was seen to have achieved the equilibrium 
state 40 element heights downstream of the element (Figure 
2-8) . At this location the boundary layer was approximately 
three times the element heights thick. Also, the sequence 
of profile shapes downstream of reattachment is just the 
opposite to that found in the case of approaching separation 
in a boundary layer flow. These experimental results were 
obtained for a strong or overwhelming perturbation, 

(6Q/h 1) , and the initial boundary layer had no influence 

on the results . 

Good and Joubert [5] experimentally obtained the 
static pressure variations across the separated region 
behind a bluff plate at various stations (Figure 2-13) . An 
important feature of these profiles is. that the locus of 
their minimums is approximately parallel to the ground and 
that it traces the region of high turbulence which initially 
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Figure 2-11. Turbulence profiles behind roughness elements [10] 









Figure 2-13. Static pressure 
a bluff plate [2]. 






coincides with the dividing mean streamline but continues to 
remain constant in height from the ground towards the end of 
the wake. The pressure losses associated with these minimums 
in the profile are believed to balance the large gradients 
in the transverse Reynolds stresses. 

III. REVIEW OF ANALYTICAL MODELS 

There is literature available regarding the 
determination of base pressure for two-dimensional steady 
(without vortex shedding) base flow. Tanner [17] has 
compared the theories of Chapman, Korst, Nash, and Kirk. 

The theory of Tanner has been discussed in comparison to 
these theories [17]. 

Most of these models assume very thin boundary 
layer at separation so that the assumption of similar 
velocity profiles is valid. In order to take the effects of 
boundary layer thickness into account, Kirk assumes that the 
shear or the mixing layer behaves as if it started some 
distance ahead of the base instead of starting at the 
separation point. It was seen that the base pressure calcu- 
lated thus will be higher than that for vanishing boundary 
layer thickness. Also, neglecting the curvature of stream- 
lines at separation (true for an "ideal" jet boundary) , the 
pressure of the uniform stream ahead of separation and the 
pressure at the outer edge of the mixing zone are the same. 

In particular, in the calculation of base pressure Chapman 
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assxames that a balance exists between mass flow drawn from 
the base region by the mixing layer and mass flow reversed 
back into the base region by the pressure rise through the 
reattachment zone, and also that the compression through the 
reattachment zone is isentropic along the dividing stream- 
line » Owing to the former assumption, for steady flow 
without bleeding into the wake region, the dividing stream- 
line, as calculated from the mixing layer theory, must also 
be a dividing streamline at reattachment. Nash improved the 
reattachment criterion by taking into account the fact that 
the pressure at the reattachment point, p^ , is not equal to 
the final recovery pressure, p^ , far downstream as had 
previously been assumed, but is actually less than pj^. 

Nash's reattachment parameter, 


N 


Pr'Ps 
Pi-Pb ' 


where p^ is base pressure, was determined experimentally for 
a back step for Mach number equal to zero to be equal to 1,6, 
and not equal to unity. 

The integral theory of Tanner predicts base 
pressure for subsonic flows and could be extended to apply 
in supersonic flow alsOo The main interest in Tanner's 
model as far as this investigation is concerned is not in 
how the base pressure is calculated, but in the way the wake 
region is modeled. Referring to Figures 2-14 and 2-15, a 
mixing process is assumed to occur in the region from the 
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Figure 2-14, Flow model for derivation of the 
second outflow function [17] , 
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R Reattachment point 
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^ from the dead- air region ceases 

(^. Outer boundary of mixing region 
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0 Particular strecunline cutting 
boundary 1 at section I 



Figure 2-15 • The subsonic steady base flow past 
a blunt body; mean streamlines and pressure 
distribution, schematically [17]. 
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separation point A to the section I which corresponds to 
constant pressure mixing between a uniform external strecim 
and fluid at rest. Due to this mixing, fluid is withdrawn 
from the dead-air region „ This phenomenon will be called 
"outflow from the dead-air region." The dividing or 
reattachment streamline separates the outer flow fluid 
from that withdrawn from the dead-air region. Experimental 
results [18] show that the velocity and hence the pressure 
at the boundary , which separates the external flow from 

the mixing region, tend to remain constant approximately up 
to a section midway between the separation and the reattach- 
ment points. At section I the mass flow between the 
boundary through which the mass outflow from the dead- 

air region takes place and the dividing streamline has 
its maximum value because then the backflow into the dead- 
air region begins. This maximum value is given by 

Zt 

m = / pUdz (2.1) 

^ z 

r 

Tanner relates m to the pressure drag and hence the base 

a 

pressure and proposes the following expression for the 
velocity profile in the mixing region (Figure 2-16) ; 

u _ 1„ . 

U2 ^ ^ ^ 2^^ ' 

where U 2 is the outer-flow velocity. This profile is an 
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Figure 2-16, Velocity profile in the mixing 

region at the end of the outflow region [17]. 
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approximation to the error function profile which is 
generally valid in a mixing layer. 

Considering the static pressure throughout the 
control voltame ABCA in Figure 2-15, page 32 (the flow 
approximates the flow in a free- jet boundary) , as constant, 
the momentum equation reduces to 

z p z 

/ y pU dz - U, / ^ pUdz == 0 . (2.2) 

z z 

r r 

Making use of Equations (2.1), (2.2), and some 
empirical parameters such as h/d. Tanner calculates what is 
called the second outflow function o Considering the region 
between I and R, Tanner then calculates the first outflow 
function. The intersection of these two functions gives the 
base pressure coefficient. However, the parameter, h/d, 
which is closely related to the spread rate parameter, o 
[19 through 25], as used in other theories, cannot be 
determined theoretically . As the spreading of the turbulent 
mixing zone is linear, a will be constant at constant Mach 
number. Different values of a have been used in the 
literature. For incompressible flow, a as a function of the 
wedge angle, 4> (for flat plate perpendicular to the air 
stream, 4> = 180®; for flow past a back step, 4> = 0°) , is 
shown plotted in Figure 2-17 (experiments conducted by 
Tanner [17]) . To account for compressibility, a is given by 
the relation, ct/Oq = 1 + 0.23 M 2 , where is the value of 


35 




Fiqure 2-17, Influence of wedge angle on the 

similarity parameter of the mixing region for 
incompressible flow [17]. 
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the spreading rate parameter corresponding to incompressible 
flow, and where M 2 is the Mach number at the outer edge of 
the mixing layer* 

The two-dimensional problem of constant pressure 
mixing of laminar two-parallel streams has been considered 
by Yen [26]. The interface velocity and the location of the 
interface in the mixing region has been determined by 
considering the conservation of z-direction momentum through 
a control volume analysis. Due to mixing, there arises a 
transverse force acting on the dividing wall. For a free 
stream velocity ratio of 0.5 and 0, the interface deflects 
towards the higher and the lower velocity streams, 
respectively . 

Kubota and Dewey [27] used a momentum integral 
technique to find the velocity profiles in the constant 
pressure laminar free shear layer. The shear layer was 
divided into two parts , one above and the other below the 
zero streamline; and separate polynomial and exponential 
expressions were used to describe the velocity profiles in 
these two regions. Closed form solutions of the momentum 
equations were obtained. 

The constant pressure laminar mixing problem was 
solved assuming the boundary layer approximations to be 
valid for the free shear layers with finite initial 
thickness. The continuity, momentum, and energy equations 
were reduced to the incompressible form by using Howarth 
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transformation. By using suitable matching conditions the 
momentum equation is decoupled from the energy equation. In 
the transformed coordinate, 

z ^ 

z = / (p/p )dz , 

0 ® 

and the streamwise coordinate, s, the equations to be solved 
are of the incompressible form. 


8u , 9v 
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where s and z are the distances measured along and normal to 
the dividing streamline, and p^ and are the values of 
density and kinematic viscosity in the undisturbed stream. 

The velocity profile is represented by an analytic 
function containing several parameters as functions of s. 

The total number of boundary conditions applied at the 
extremities of the shear layer, at z = 0, and the moment 
equations obtained by multiplying the streamwise momentum 
equation by u^ (j - 0,1,2,,..,) must equal the number of 
parameters appearing in the velocity profile. As the down- 
stream velocity profile is quite different from the initial 
profile due to the non-similar growth of the free shear 
layers with a finite initial thickness, it is difficult to 
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represent the profile as a whole by a single expression. 
Hence the layer is divided into two regions, one above the 
zero streamline and the other below it. The profiles for 
these two regions are assumed to be of the form, 


and 
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respectively; < 5 ^^ and 62 are the upper and lower shear layer 
thicknesses, respectively; m and n are so chosen that the 
profile is realistically approximated. The profile 
parameters, Aq , a^^ , Bq , bj^ / / and 62 are determined 

by the boundary conditions at the lower and upper 
extremities of the shear layer, the continuity of the 
velocity and its derivatives at z = 0, and the moment 
equations. The possible boundary conditions for the outer 
and inner boundaries are, respectively, 


a = 0 n = 1 

} f =: 1 , f* = f" = •••• = 0 

0 = 1 ri “ 

and 


a = 0 h = -1 

}g=0,g*=g"=***«=0, 

0=1 h 
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where the primes denote derivatives of f and g with respect 
to n and h, respectively. The matching conditions at z = 0 
are given by 




z-O 


(^) (—1?) ; j = 0,1,2,* 

2 fih-' 


For a quadratic profile, a = 0 and m = n = 2; 

letting j = (0,1), it was shown that (u/u ) was 0.596, 

e z~ u 

which is very close to the exact value of 0,587 found by 
Chapman [28], Using an exponential profile, (a ~ 1) , 
letting m = n = 0, and j = (0,1), (u/u ) __« was found to be 
0.618. 

At the separation point, s = 0, z=0 ” 

^2 ~ equal to the initial shear layer thickness, 

6q, In the quadratic and the exponential profiles only one 
parameter, (z/6q) , specifies the initial profile. If twor or 
more such parameters were available, the effect of the 
initial profile on the rate of growth of (u) would be 
known. Figure 2-18 shows the effect of the shear- layer 
profile on velocity along the zero streamline. 

The Karman-Pohlhausen method for two-dimensional 
laminar flows does not give reasonable results in the region 
of adverse pressure gradient, and particularly when 
separation occurs^ Between the separation and reattachment 
points this method fails to represent any real situation. 
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Figure 2-18. Effect of shear layer profile on velocity 
along zero streamline [27]. 



Where the pressure gradient is zero, the Karman-Pohlhausen 
method must yield an attached Blasius-type velocity profile. 

In Figure 2-19 reattachment will be predicted upstreeun of B, 
the region of zero pressure gradient, where in actual fact 
reattachment occurs downstream of B. An integral method 
first proposed by Walz [29] and modified by Tani [30] has 
been demonstrated by Lees and Reeves [31] to produce 
profiles with reverse flow even for zero pressure gradient 
and, therefore, could predict the behavior of separated flows. 
Also, this method does not require the empirical data which 
other methods require, for example, modified Crocco-Lees 
theory [32 ] . 

Modeling the flow separating from the sharp edge 
of a back step as a free- jet boundary between the outer flow 
with initial velocity, u^ , and a lower stream with a 
velocity of zero, and assuming zero initial boundary layer 
thickness, similar velocity profiles are given by the 
expression for the velocity profile in the region of a jet 
interaction [22] as 


u = (1 + erf C) t 

where the transverse coordinate. 
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Figure 2-19. Typical static pressure variation 
on a flat plate in supersonic flow With 
separation induced by an incident shock wave. 
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and 


5 2 

erf 5 = — / e"^ dt 



The spreading parameter, a, is an empirical 
constant. Various values of a ranging from 9 to 15 have 
been used for back steps for incompressible flow. Plate 
[23] used a = 14.5; and the x-axis was taken as the locus of 
points where u/u^ = 0,5, and the z-axis as perpendicular to 
X, Comparison of Plate's solutions with the experimental 
data (Figure 2-20) shows some disagreement at the negative 
value of ^ where the lower boundary begins to influence the 
data. This is attributed to the theoretical model which 
assxames that the lower velocity stream extends to minus 
infinity. The theory of free-jet boundary yields a Gaussian 
shear stress distribution, which is in good agreement with 
the experimentally measured shear stress distribution [33“] . 

To make use of the error function profile, the 
curve given by u/u^ = 0.5 must be known either experimentally 
or theoretically, A theoretical prediction (work in progress 
by Bitte and Frost at The University of Tennessee Space 
Institute) can be made on the basis of free streamline 
theory, which is based on the velocity distribution in the 
undisturbed boundary layer, the shape of the obstacle, and 
some empirical input related to pressures on the obstacle 
and the reattachment length. 
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Figure 2-20. Correlation of velocity in the shear layer with 
error function profile [2]. 



Camarata [21] correlates the velocity profiles 
obtained from experimental data in the similarity region 
close to separation with the power law profiles. 


u 
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with n varying from 6 to 9. The correlation parameter is 


If (U/Ug = 0.5) , 

where 

(f> = u/u^ 

and 


C = z/6q 

u^: free stream velocity, 

initial boundary layer thickness. 

With this correlation and the two- layer (outer and 
inner mixing zones) model for the velocity profile in the 
shear layer, an integral procedure [21] which does not need 
any specification of the shear stress distributions is 
readily applied to calculate the profiles. 

The pre~asymptotic region of the free shear layer 
is divided into an upstream segment and a downstream segment 
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(Figure 2-21) . The profile in the shear layer in the down- 
stream segment is an error function profile, 




In the upstream segment, the shear layer is partly a 
truncated error function profile for the inner mixing zone 
and partly a power law profile for the outer mixing zone, 
the free-stream edge of which does not appreciably change in 
height above the separation point at any x- location. The 
reason for this is that the initial boundary layer has a 
very slow growth relative to the rate of spreading of the 
mixing layer. Thus ~ constant for all x stations. Since 
the whole initial boundary layer cannot be affected 
immediately by the abrupt removal of the wall, it is 
physically realistic to subdivide the shear layer as 
mentioned above. 

Referring to Figures 2-22 and 2-23, an empirical 
correlation curve specifies the slope of the velocity 
profile in the pre- asymptotic region at the point, Az, where 
u/u^ = 0.5. Therefore, 

1-n 

|H (u/u^ = 0.5) = (A.) " (2.3) 

The left-hand side of this equation is known; and for a 
given x, n is known. Hence Equation (2.3) will yield Az. 
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Figure 2-21, Flow model for pre- asymptotic iregion [21], 


z 



Figure 2-22. Pre- asymptotic velocity profile 
description [21]. 
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Figure 2-23. Variation of velocity profile slope at 
the half-f ree-stream velocity point with 
distance from separation [21] . 
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aAz 
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(2.4) 


Also, 


du u 

erf / A _ e a 
—j- '■ (Az) — — — 


du 

dz 


(U/Ug 


0.5) . 


For a given x, a can be iteratively calculated from 
Equation (2.4). Therefore, the error function profile is 
known . 

Conservation of momentum equation yields 

j 2/n^ u^ « ^ 

/ u^(z/6q) ^ dz = / -^(1 + erf ^) dz , (2.5) 

0 —00 


where n^ (the value at separation) has been taken as 7 [21] . 

The preceding equation can be solved for z^ 
iteratively. Hence the velocity profile is completely known. 
Comparison of experimental data with the solutions is shown 
in Figure 2-24. 

Far downstream from the asymptotic region, which 
is not a similarity region, no correlation is available; and 
hence a has to be specified. From the theoretical estimates 
of earlier authors [21] o in this region is approximated by 
the relations, a = 12 + 0,8 free-stream 

Mach number. Hence the transition between the asymptotic 
and the pre-asymptotic regions occurs at the point where the 
values of a obtained from the two criteria are equal. 
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n =s 7 

6q = 2,7 in. 



4 ) 

Velocity ratio, u/u^ 


Figure 2-24. Comparison of theory with incompressible 
measurements of Charwat and Der [21] • 
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Korst and Chow [19] have compared various 
analytical and experimental constant-pressure free shear 
layer similarity profiles by the use of spread rate 
parameters. It was found that in addition to the flow 
conditions, a also depends on the analytical model, the way 
matching is done with the experimental data, and the 
viscosity model of turbulence employed. Although matching 
along the dividing streamline could be done, matching at the 
half -free-stream-velocity point, u/u^ = 0,5, was observed to 
be more realistic, as shown by Figures 2-25 and 2-26. 

Various profiles agree very closely at this point. 

The initial conditions like the initial boundary 
layer thickness, the initial velocity profile, and the 
initial distribution of turbulence intensity influence the 
flow downstream in the free shear layer. As long as the 
change in the pressure gradient is small, the turbulent flow 
equations can be laminarized to some accuracy through the 
use of an eddy transport coefficient concept. Turbulence 
shear stress, in analogy to laminar shear stress, is written 
as 

9u 

^ - 3 ^ ' 

where e is the eddy viscosity. Various eddy viscosity 
models have been examined [20] . These viscosity models as 
applicable to incompressible flow are as follows: 
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0 0.2 0. 4 0.6 0.8 1.0 


Figure 2-25. Comparison of different 
velocity profiles [19] . 
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- kj_0U 


= k2(x + Xq)u^ 


s k u^/ ( ^^ ) 


* k4(x + Xq) 


/ 


3u 


(2.7) 

( 2 . 8 ) 

(2.9) 


where ^ 3 ^ f ' ^3 ' ^4 empirical constants. The 

constant x^ is the shift of the virtual origin due to the 
presence of an initial boundary layer. It accounts for the 
initial thickness. Subscript D refers to the dividing 
streamline , and e indicates the outer edge of the mixing 
layer (Figure 2-27) . 

In Equation (2.9) Prandtl's mixing length, 
expressed by v1^ (x + ^ is assumed proportional to the 

streamwise distance, x. Equations (2,6) and (2.7) are 
similar to Prandtl's simplified hypothesis [22] in which e 
is assumed proportional to the width of the mixing layer. 

Nash [24] used Equation (2.8) in his study of the pre- 
asymptotic shear layer. 

At the separation point a singularity exists, and 

there is a small neighborhood where a discontinuity in shear 

stress occurs. This neighborhood for flat plate has been 

-3/4 

found to be 0(LR ' ), where L is the length of the plate and 
R is the Reynolds number based on L. Neglecting this 
extremely small region, laminarized boundary layer equations 
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InvlBcld flow 



Figure 2 - 21 . Free shear layer model [20] • 



describe sufficiently accurately the similar and nonsimilar 
flow regions in a free shear layer. The effect of the sharp 
turn of the attached boundary layer flow on the initial 
velocity profile is neglected. Therefore, the attached 
boundary layer profiles are the initial conditions for the 
separated shear layer. 

Values of the empirical constants evaluated after 
comparing the experimental data with the calculated velocity 
profiles 120] were found to be as follows; = 0.6; 

■ 0.002; k^ = 0,015; k^ = 0.0575; x^ = 30 6 q , where 6^ 
is the momentum thickness at the shifted origin. The 
similarity parameter, a, was taken as 9.42. 

The model as described in Equation (2.9) does not 
give results in agreement with the other three. This is 
because it yields a different viscosity distribution. 
Viscosity increases from zero to a maximum inside the shear 
layer eind then decreases to a value of zero at the upper 
edge of the mixing layer. But in the other three models the 
viscosity is constant across the layer. For the sake of 
simplicity these three models are well suited for calculating 
the meem quantities in turbulent free shear layers. However, 
although the empirically determined values of the k's give 
good correlation in one region, they may fail in other 
regions of the same flow field. Hence they are not universal 
constants for the whole flow region. 
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In the case of mixing in a free- jet boundary^ 

streamlines have a curvature which becomes important and 

should be taken into account. Uchida emd WataneQ^e [34] 

defined streamline coordinates in a curved flow field by 

d s 

a and 3 and their extension parameters by h^ « and 

Jin 

h^ * *33 ' s ^ represent arc lengths along and 

normal to the streamline/ respectively (Figure 2-28). - 

h^ and h^ satisfy Gauss* orthogonality equation/ 


3 , 1 ^*^6 


('T 

a 


3,1 '»»», 


l\ + _ n 

3a'h 3a ^ 3B^h« 33 ^ ' 


with 3 = constant chosen as streamlines; 


^ 

hg 33 hg 


-f H = o , 

a 


R = 


u i 
s s 


The concentration of the mass flow in the mixing 
zone expressed by 3 is lower than that expressed by 4* the 
Cartesian coordinates , and hence 


3u 

33 3z 
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Figure 2-28. Streeunline coordinates system [34]. 
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Applying boundary layer approximations to the Navier-Stokes 

equations in the streamline coordinates, two important facts 

become apparent [34], Firstly, the nondimens ional pressure 

gradient, , counterbalances the centrifugal force; and 

secondly, the effect of curvature becomes important when it 

1/2 

is of the order of R ' . In making the boundary layer 

approximations , orders of magnitude for various quantities 
are assumed to be 


0(S) = 

0(h^ , h^ , u , p , a) = 1 . 


At 3 the solution is matched with the outer flow of 

velocity, u , Defining the variable [35], 
e 


X 


“a , 1 . 

vT' ~ 


$ 


1/h^ becomes equal to u^ on the outer boundary where X = 1, 
Thus the boundary conditions can be expressed as 


at 3 


= 1 and = 


t 


at 3 ; X = 0 . 

In deriving similar solutions, a similarity 
variable is defined as 
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n = /rToF , 


which takes into account the fact that the width of the 
laminar mixing zone increases as corresponding to x^'^^ 

in the usual case. 

Assuming no pressure gradient along the zero 
streamline for incompressible laminar mixing in free-jet 
boundary, the authors theoretically generate shear stress 
along = 0 as a function of the curvature parameter, C, 

along = 0 (Figure 2-29) , where 


_ 1/2 




H is a function of integration. 
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4 -0. 3 -0.2 -0.1 0 0.1 0.2 0.3 0. 

C ^ 

Figure 2-29. Shearing stress along zero 
streamline [341. 




CHAPTER III 


SOLUTION OF GENERAL NAVIER-STOKES EQUATIONS AS 
APPLICABLE TO ATMOSPHERIC FLOW 


Governing Equations for the Turbulent Atmospheric Flow 

Equations of motion for a steady, two-dimensional, 
incompressible Newtonian flow in the Cartesian coordinates, 

X and z (Figure 3-1) , are 


9u _9v 
3x 9 z 


0 


(3.1) 



9u _ 1 9p 

'^•5^ " ■ p 3$ 


+ 


i -^[u (2—) 1 

p ax^^eff^ 9x^^ 


+ 


9 j. , 9u 
9z^^eff ^9z 


+ 



(3.2) 


9 


+ 



1 9p 

p 


+ 


1 9 , /o SVv T 

p 9z ^^eff ^^9z^ ^ 


9 p p9u 9v 

dx^^eff^Jz 9x) ] , (3.3) 


where u and v are the ensemble averaged components of 
velocity along the x and z directions. 

In this study these equations are assumed to 
govern the neutral atmospheric boundary layer flow. The 
Coriolis forces induced by the rotation of the earth [36] 
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are considered negligible in the lower regions of the 
atmosphere to which the present problem is confined. 

Boundary Conditions for the Atmospheric Shear Flow 

The thickness of the atmospheric boundary layer is 
of the order of 1,000 meters. Thus the outer boundary 
condition for the problem under consideration is specified 
from within this shear layer. At the ground no-slip 
condition must be satisfied. In an equilibrium turbulent 
boundary layer the logarithmic velocity profile is given 
by [36] 


u 



z+z, 


«,n 


(3.4) 


The introduction of the surface roughness parameter, z^ , 

ensures that the mixing length does not vanish at the 

surface; and hence the turbulence is not zero at the surface 

* 

(z = 0) . The friction velocity, u , is taken to be constant 

from the assumption that the boundary layer is a constant 

stress layer = t ) . 

■* xz w 

Eddy Viscosity Model 

For a laminar incompressible flow molecular 
viscosity, y, can be assiimed constant; but in a turbulent 
flow the effective viscosity, / comprised of turbulent 

and laminar parts, is given by 

‘'eff = ^ 
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where , also called the eddy viscosity/ is a function of 
the velocity gradients. From Prandtl*s mixing- length 
hypothesis , 


y 


t 



3u 


(3.6) 


with the mixing length. 


£ = k(z + Zq) . (3.7) 

Equation (3.7) is valid in the immediate vicinity 
of the wall. Away from the wall the mixing length distri- 
bution depends on the flow situation. Equation (3.6), 

applicable to boundary layer type of flow, gives a zero 
3 \i 

viscosity for = 0 which is not necessarily true for 
recirculating flows where the velocity gradients, and , 
may be comparable in magnitude. Hence the eddy viscosity 
model used in the present study is taken as 

2 2 

Ut = + (|^) } (3.8) 

Essentially, the turbulent flow has been treated as one with 
turbulent viscosity as its transport property which is a 
function of the location of the flow region. For recircu- 
lating flows, however, where the relations between stresses 
and velocity gradients are complicated and not yet known 
[37], the eddy viscosity model is unsatisfactory. Therefore, 
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one is lead to a more detailed model of turbulence developed 
by Prandtl and Kolmogorov, 

Prandtl-Kolmogorov Model of Turbulence 

According to the Prandtl- Kolmogorov model, here- 
after referred to as the TKL model, the turbulence kinetic 
energy, k, and turbulence length scale, I, characterize the 
state of local turbulence. This model provides a better 
formula for determining the turbulent viscosity since for 
recirculating flows where the mean velocity gradient 
vanishes, the turbulence energy is not necessarily zero. 

Algebraic equations for the transport properties . 
The TKL model relates k and I to the effective viscosity by 
the relationship. 


Ueff = % O' ' 

where k and £ are the local values , and is a function of 
turbulence Reynolds number which is defined as 


R 


t 



(3.10) 


where the Prandtl-Kolmogorov formula for turbulent viscosity 
[37] is given by 


y 


t 





(3.11) 


it can be argued that when R^ is small, turbulence is 
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negligible and the function, , tends to 1/R^* On the 
other hand, when is large, it takes on a constant 
asymptotic value, C . Researchers have used different 

^oo 

asymptotic values for C for different flow situations, but 

r* 

no "universality" in its value has been found so far. The 
knowledge of the behavior of as a function of R^ in its 
intermediate range is also not satisfactory [38]. 

From Equations (3.5), (3.9), (3.10), and (3.11), 

the following relation emerges: 

= 1 + 1/R^ , (3.12) 

which implies that C =1. This value has been used in 

^co 

different flow situations [38] , but in recirculating flows 

it is different from 1. Further discussion on C will 

^00 

appear in the sections on selection of constants and their 
parametric study. In general, 

C = C + 1/R . (3.12a) 


Equations to be Solved 

Differential equations for the vorticity, w, and 
the stream function, ip , are obtained from Equations (3.1), 
(3.2), and (3.3) by making use of the following transfor- 
mations . 
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y-axis. 


Define a two-dimensional vorticity about the 


, _ 8v 3u 


(3.13) 


and a stream function, \}j, such that 


u 


1 

Ip Tz ' 


(3.14) 


V 


1 dip 

p JSi • 


(3.15) 


Introducing ip into Equation (3.13), eliminating 
pressure from Equations (3.2) and (3.3), and introducing o) 
and ip into the resultant equation gives the Poisson equation 
for the stream function and the transport equation for 
vorticity, respectively. Both of these equations can be 
expressed in the following general form of the elliptic 
equation [38] , 


r r 9 / j,dlp 


)}] 


- + <j = 0 , (3.16) 


where the coefficient functions, a, b, c, and d, and the 
dependent variable, (p, are given in Tables 3-1 and 3-2. In 
Equation (3.16) the terms in the first pair of brackets are 
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TABLE 3-2 


EXPRESSIONS FOR THE SOURCE TERI4, -d 
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the convection terms, those in the second one are the 
diffusion terms, and the last one is the source term. 

Gosman, et al . [38] , derived a differential 
equation for k, and Rotta [39] proposed a differential 
equation for both of which can be expressed in the form 
of Equation (3.16). 

Referring to Tables 3-1 and 3-2, the term in 
brackets multiplying in the source term for the 
k-equation is the production term; the second term repre- 
sents energy dissipation. The term containing C in the 

s 

Jl-equation represents the rate of increase or "stretching" 

of the length scale, and the second one represents the 

decrease or "breaking" of the length scale. The coefficients 

C , and C , are believed to behave with R. like C , given in 
a s t y 

Equation (3.12a) [38]. 

Solutions for laminar flow were determined at 
different Reynolds numbers solving the two partial differ- 
ential equations (pde) for w and Turbulent flows were 

solved using the eddy viscosity model introduced into these 
equations for intermediate values of Reynolds number. For 
higher Reynolds-number turbulent flow, the TKL model was 
employed which introduces the two additional partial differ- 
ential equations for k and Z in the analysis . 

Derivation of the Finite-Difference Equation 

The finite-difference equation (fde) has been 
derived by integrating the general pde given by Equation 
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(3.16) over a finite area, assuming that the distribution of 
the variables between the nodes of the grid is known. 

Integration of Equation (3.16) over the rectcuigular 
domain (sides of the rectangle lie midway between the 
neighboring grid lines) shown by dotted lines in Figure 3-1, 
page 65, will give 

X z 
w s 

X z 

X z 
w s 

X z 

+/®/^ddxdz=0 

X z 
w s 

To illustrate the principles involved in inte- 
grating the convection terms , one of the four integrals , 
namely , 


s e 

is integrated as follows : 

Assuming that and 4 / are well-behaved functions 
and that an average value of (J)^ exists, which can be given 
by 
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J = 


z 

z 

s e 

(||) dz 


the integral takes the finite difference form 


a?g(0-„g - 'l'se> • 


Assuming that (p within each rectangle in Figure 
3-1, page 65, has a constant value equivalent to that at the 
centrally enclosed node, P, and that the takes on the 
4>-value possessed by the fluid upstream of the e-face of the 
rectangle, the upwind differencing effect is introduced into 
the integral by the following finite difference form: 




(ll; )- ]p 

^ne ^se ^ne se 


1 . A / ^^ne^’^se^ I ^ne^^se I -L 1 

■t + K n 


If the flow direction is from P to E, then \Jj -ip is 

ne se 

positive and hence 



If the flow direction is from E to P, then 
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Also, assximing that the ip-value at a particular 
corner of the rectangle is equal to the average of the four 
neighboring node values, , for example, can be given by 




+ 1'P + ’^E 


ne 


The remaining convection terms are integrated similarly. 

Considering the first term in the set of four 
diffusion terms, namely, 


z 


it is assumed that b and (ccj)) can be approximated as varying 
linearly with x so that 


and 




d 

JSF 


(C(^>) 


^E ■ ^P 


Thus the finite difference form of the first diffusion term 
is 


^E ^P ^E^E ” ^P^P 

5 Xj. - Xp 


(z - z ) 
n s 
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with the other diffusion terms having similar expressions. 
Finally, the last term, 

X z 
e n 

/ / d dx dz , 

X z 
w s 

is integrated by assuming that the value of d is constant 
over the area of integration and is equal to that at the 
centrally enclosed node, P; thus the source term becomes 


d_ (x - X ) (z - z ) . 

P e w n s 

On rearrangement and simplification the following 
fde is obtained from integration of Equation (3.16) [38]: 

Aj,(*p - (fig) + - ♦„) + A^(*p - + Ag(<t>p - .j)g) 

®E^°e'*’e ” ■■ °p‘*’p^ " ■ °p'*'p* 

- Bs(°s+S - °P*P> + Vp = ° 


where the A*s, the coefficients in the convective terms, are 
given by 




'^se'*''^s‘''^ne"'^n 




/ 
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Ag = ^'I'sw''’’^w“’^se”'^e^''’ '^sw''’’^w”’*’se”’^e 


] , 

] , (3.17a) 


and the B*s, the coefficients in the diffusion terms, are 
given by 


and 






W 




N 


B... - 


bg+bp 

^n’^s 

— 3 — 

Xg-Xp 


^n”®s 

—J— 

Xp-Xj, ' 


^e'^w 

— 3™ 

^n"^p 

bg + bp 


4 



Vp - (— ^) (— ^ 


(3.17b) 


(3.17c) 

3(j) 


The first order derivatives, for example, of 
the dependent variable, (j), are approximated as follows: 




^N"*^S 


* <V^s^^g 

dz 




(3.17d) 


It may be noted that in deriving the coefficients 
in the convective terms , use was made of "upwind 
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differencing" for first order terms so that the A*s are 
always positive. This implies that the transport! ve 
property is advected in the direction of velocity which is 
realistic. The stability of the solution procedure improves 
considerably using upwind differencing. The coefficients in 
the diffusion terms, the B*s, are dependent on the grid 
spacing and are also all positive. One of the assumptions 
made in reducing the diffusion terms to their finite 
difference form is that b and (c(|)) vary linearly with x over 
the domain of integration (Figure 3-1, page 65) , 


Complete Successive Substitution Formula 

Equation (3.17) can be written as 


where 




Cg = (Ag + BgCg)/Z: AB , 


‘'N ” ^ ' 


Cg = (Ag + BgCg)/2: AB , 


D = - dpVp/Z AB , 


and 


Z AB = 




where the C*s are all positive. 
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This sxibstitution formula can be applied at every 
interior node in the flow field. The simultaneous nonlinear 
equations for oi, ip, k, and Z given by the general equation 
(Equation (3.18)) are solved by a Gauss-Seidel iterative 
procedure . 

Equation (3.18) for k can be written as 
where the source term, S, „ , is given in Table 3-2, page 72. 

K. fir 

This formula has been modified [38] by rearrangement to give 
(refer to Table 3-2) 


, ^E^E ^skt,P^p/^ ^ 

kp = 7“ — ' 

^ ''' Iab /^)p 


. (3.18a) 


For Z, Equation (3,18) becomes 


^P ^E^E 


+ (pk^/2 Cg - I jIp , 


which can be rearranged to yield 


S = 


" ^ rlB<\ktV^>P 


(3.18b) 
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The advantages in using Equations (3.18a), 

(3.18b), and the implicit substitution formula for o) 
(Equation (3.31)), which is given in the section on boundary 
conditions, will be discussed in the section on convergence 
of the solution procedure. 

Setting Up and Discussion of the Boundary Conditions 

Referring to Figure 3-2 and considering the flow 
field enclosed in a control volume 1-2-3-4-5-6-1 with a unit 
depth perpendicular to the plane of the paper, flow enters 
at 1-6, the inlet, and passes over the back-step 1-2-3 and 
leaves at 4-5, the outlet. The lower boundary 1-2- 3- 4 
comprises the top face of the step 1-2, the base 2-3, and 
the ground 3-4. Surface 5-6 is referred to as the upper 
boundary. The Cartesian coordinate system is shown in 
Figure 3-2 with the origin at station 3. 

The boundary conditions (b.c.) prescribed on <p are 
either of flux-type or of normal-gradient- type, except at 
the inlet and the outlet where the (p values prescribed are 
also related to the values of the horizontal component of 
velocity for all locations on the boundary. The normal- 
gradient- type boundary conditions retard the rate of 
convergence in comparison with the flux- type boundary 
conditions . 
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Algebraic equations for and substitution formulae 
derived from the boundary conditions « At the inlet and 
outlet, logarithmic velocity profile (Equation (3.4)) 
typical to an equilibrium atmospheric boundary layer is 
prescribed with different shear stress velocities. At the 

outlet, the reference velocity, u^ , is prescribed at the 

* 

right-hand corner of the upper boundary so that u^ is given 
by 


* 



Zn 




(3.19) 


Considering the upper boundary to be in the undisturbed flow 
and applying the continuity equation, 


^0 * 
u. 

/ —in 

K 


Z+Z, 


dz = 


“O * 

f ^ Zn ^ dz , 

K Z, 


0 


0 


so that the friction velocity at the outlet becomes 


* 

u. 

1 


.11 


0 


0 

h«-h+z, 


ho”h+Zo 

[Ji )i{Zn — 

^0 ^0 


(3.20) 


- 1 ) + 1 ] 


Equations (3.19) and (3.20) imply that iff is constant along 
5-6. 
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Boundary conditions for the stream function . The 
description of the boundary conditions on ip that were used 
in the numerical solution of the problem is as follows: 

The boundary conditions on different boundaries 

are : 


Inlet 




S’l' 

h ^ 


pUiZQ Z-h+ZQ 


z-h+z 


(in 


0 


- 1 ) + 1 ] 


( 3 . 21 ) 


Outlet 


1 /) - / 
0 


dZ 


PUqZq Z+2 

dz = -_[-_(in 


z+z, 


- 1 )' + 1 ] 


( 3 . 22 ) 


Lower Boundary 


ip = 0 (no-slip condition) . 


Upper Boundary 


^ = \pQ - constant . 

Unless the outlet is sufficiently far downstream, 
fixing ip there through the logarithmic velocity profile may 
cause the solution to be rather unrealistic. In the present 
problem the outlet has been taken far enough downstream that 
the logarithmic velocity profile is a reasonable approxi- 
mation to the real flow situation. In order that the outlet 
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b.c. could exactly represent the true physical situation 
which theoretically occurs at infinity, use would have to be 
made of a transformation like n - which maps the semi- 
infinite axis, 0 ^ X < 00 , into 0 < n < 1 [40], As the 
asymptotic behavior of the outlet velocity profile is that 
of the logarithmic velocity profile, this transformation 
could hold promise; but care would have to be taken in 
choosing the proper grid distribution, especially in the 
recirculation region. Using this mapping also introduces 
some more gradient terms in the source term of Equation 
(3.16), which can influence the behavior of the iterative 
process as far as convergence is concerned. 

Roache and Mueller [411 used a floating b.c. on 
at the outlet; 

A = 0 

9x 

with (3.23) 


which amounts to extrapolating ip linearly as 

The preceding relation is valid for a uniform grid near the 
outlet in the x-direction. The authors [41] report that 
Equation (3.23) failed for low Reynolds numbers as it gives 


an abrupt variation of o) near the outlet. However, for high 
Reynolds numbers it gives good results. To use these 
boundary conditions at the outlet, u^ is taken as the 
reference velocity at the inlet at the upper boundary so 
that 


* 

u. 

1 


h^-h+^ 


in 


(3.24) 


This specifies the inflow b.c. on ip given by Equation (3.21), 
The b.c. for the outlet given by Equation (3.23) along with 
Equations (3o24) and (3,21) have also been investigated for 
the problem under consideration. Further discussion of 
these boundary conditions appears in the section on 
discussion of results. 


Boundary conditions for the vorticity . The 
vorticity boundary conditions as prescribed on the various 
boundaries are given belowo 

Inlet 

At the inlet the following equation for (c was used 
with Equation (3.21) for with either of the equations, 
Equation (3^20) or Equation (3o24): 


8v 9u 
9x " 9z ' 
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or 


X3-X2 “ X2-X, u. 

rrx^-x'l ■* " K ( z+^T • 


(3,25) 


The expression for — in the above equation 

follows from the second order Taylor expansion neglecting 

.3. 

— !£ and higher order terms, whence 



(3.26) 


Equation (3,25) allows -1^ to develop as a part of the 
solution o The flow is thus not completely specified at the 
inlet lest the elliptic nature of the problem be restricted 
[42] , 


Outlet 


With the b.c. on ijj given by Equation (3.22), the 
following b.c. on (ii was employed: 


^IN~^INM _ ^INM~^INM-1 


^0 

tc (Z+Zq) 


(3.27) 


Upper Boundary 


_ 9v _ 9u 
^ 9x 9z • 


As 


9^ij^ 

9x^ 


at the upper boundary is zero. 


it follows that 


9z 
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Extrapolating u linearly to the upper boundary requires 


9^u n • 9^^ A 

— y = 0 , i.e., TT— = 0 


which gives the relation, 


^IN ‘^INM ' 


(3.28) 


as used by Roache and Mueller [41] • 


Lower Boundary 

Using the no-slip condition, the continuity 
equation, and the fact that v = 0 along the lower boundary, 
and by retaining the third order terms in the Taylor 
expansion for (see Figure 3-2, page 82), Woods [43] 

arrived at the second order form for vorticity at the solid 
wall. 


^NP 2 

0) = - + -^] + 0(Az"^) , (3.29) 

® p(Az)'' 

which assumes that 03 varies linearly with the normal 
distance from the wallo Gosman, etal. [38], have achieved 
more accurate predictions with Equation (3.29) than with the 
first order explicit form, 

2\l) 

oa + 0(An) , (3.30) 

® p(An)"^ 
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where n is the normal distance from the wall, and o) is, the 

- , , ' • -L S . ,,r !. 

value of vorticity at the solid wall. These two forms are 
valid for a uniform grid near the wall. Referring to 
Figure 3-2, page 82, Equation (3.29) can be used implicitly 
in the substitution foannula (Equation (3.18)); and the 
modified substitution formula for vorticity for the corner 
node NP, for example, can be written as 


^E^E 


‘^NP 


^ ^ 3,VnP 

“ p ^7777 
(Ax) 


^sV 


] + D 


1 + 




(3.31) 


If the explicit boundary condition given by Equation (3.30) 
is used, the appropriate boundary condition at the sharp 
corner 2 is not well defined, and a number of foraulae have 
been proposed [22]. Two of these are 

1. 0i2 = 0 n _ 

2, Discontinuous values of vorticity, 

-2i(jj^p/p (Az) ^ and -2if;j^p/p (Ax) ^ are used, 
depending on whether the vorticity at the node 
immediately above , 2 or at the. pne imniediately 
downstream of 2 is being calculated, 
respectively. ^ 

With the implicit scheme (Equation (3.2Q)), it is 
unnecessary to specify any b.c. at 2 because the discon- 
tinuous nature of vorticity at 2 has already entered the 
substitution formula (Equation (3.31)). Comparisons of 
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these explicit and implicit forms were made in the beginning. 
No appreciable difference was observed although the latter 
seemed to yield more realistic results. 

Boundary conditions for the turbulence kinetic 
energy . Deriving boundary conditions on k necessitates 
finding a simple expression for it. Differentiating 
Equation (3.4) with respect to z will give 


* 

3u _ u 
3z K {z+ZqJ * 


(3.32) 


Referring to Figure 3-2, page 82, the wall velocity gradient 
near the ground and the top face of the step will be and 
that near the base will be • Restricting attention for a 
moment to the regions where >> , Equation (3.6) and an 

expression originally suggested by Boussinesq in analogy to 
the expression for laminar shear stress, namely. 


T 


= y 


3u 

t 


f 


give the following relation: 


Also, 




3u 3u 
3z Iz • 


(3.33) 


eff t 
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Using Equation (3.32), laminar shear stress will 
be given by 


^ K (z+2q) ' 


Therefore , 


2 , * 

T = PU r X ' , 

eff K (Z+Zrt ) ^ K ( z+ Zrt ) 


pu u ^ .1 

( z ? 2 q ^ ^ * 


(3.33a) 


By definition. 


= pu 


Assuming that the magnitude of the shear stress is constant 
near the walls and throughout the undisturbed layer# 


eff 


eff 


w 


pu 


(3.33b) 


Equating Equations (3.33a) and (3.33b) gives 


I = [k (Z+Zq) (k (Z+2q) ^}] 

u 


1/2 


Since k(z+Zq) >> , it follows 

u 


A ~ k(z+Zq) . 


Also, 


(3.34) 


eff 


du 

= P^ff ^ • 


(3.35) 
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From Equation (3^9), 


V = c 

eff ^ • 


Using Equations (3,32) and (3.34), Equation (3.35) becomes 


Tgff = pc^k^/^u* . (3.35a) 

Equating Equations (3.33b) and (3.35a), 

.k = (u*/C . (3.36) 


Near the base, >> ; it is assumed that the 

velocity variation in v with x is also governed by the 
logarithmic velocity profile, 


V 


* 

V 

K 


x+z, 


Zn 


so that 


(3.37) 


* 

3v _ V 
'Sx K (x+Zq) 


(3.38) 


Writing the analog of Equation (3.33), 


T 


xz 



3v 3v 
3x 3x 


t 
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and following the same reasoning as before, it can be seen 
that near the base 


I K(X+Zq) , 


(3.39) 


and 


k 


* 2 
(— ) 
% 


(3.40) 


Equations (3.10), (3.11), (3.12a), (3.34), and (3.36) give 


k - [{u* 


y 

PK (Z+Z 




(3.41) 


Similarly, Equations (3.10), (3.11), (3.12a), (3.39), and 
(3.40) yield 




(3.42) 


9 H ^ 

For the inlet and outlet, >> ; therefore, at the inlet. 


■ p(z+z„) ^ C„ ^ ' 


(3.43) 


and at the outlet. 


- pTzfe^> • 

^00 


(3.44) 
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with these boundary conditions, k is fixed at the 
inlet and the outlet. But, corresponding to Equation (3.23), 
a floating b.c.. 


3^k 


= 0 




was used at the outlet, which gives 


^IN " ^^INM " ^INM-1 • 


(3.44a) 


At the upper boundary the b.c. is = 0, which in finite 
difference form is 


^JN ^JNM * 


(3.45) 


At the lower boundary the b.c. on k can be derived as 
follows. From Equations (3.32) and (3.38), respectively. 


and 


* 


u 


< (z+Zq) 

p 



* 

V 


K (X+Zq) 



/ 


where the following expressions for 
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and 



at the respective walls, obtained by using Taylor series in 

which third and higher order terms are. neglected, are 
2 2 

2i/^j^p/(Az) and 2ipj^p/ (Ax) , respectively (see Figure 3-2, 
page 82) . Therefore, the preceding two equations at the 
respective walls become 


and 


* 

u 


2kz^i1j 

p (Az) ^ 


(3.46) 


* 

V 


2kz^ 
p (Ax) ^ 


(3.47) 


Referring to the corner node NP (Figure 3-2) , Equations 
(3,41) and (3.46) give 


Similarly , 


k - r_L_{!!!o!s£ - 

s ■ 'pc (117^ *="0 




- r-i-{ 2 

PCy (AX) 2 ^^0 


(3.48) 


(3.49) 
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At the sharp corner 2, discontinuous values of k, 


^ 2K(h+ZQ)\p 


NP 


pC 


y. 


M 1 1 

(Az)^ ' 


and the one given by Equation (3.49) are used depending on 
whether the turbulence kinetic energy, k, at the node 
immediately above 2 or at the one immediately downstream of 
2 is being calculated, respectively. 


Boundary conditions on Referring to Figure 3-2, 

page 82, the distribution of the mixing length is shown 
which is based on the normal distance from the nearest wall. 
The lines 22' and 33' are inclined at 45® from the positive 
x-axis. This distribution fixes the mixing length (1) on 
all the boundaries for the TKL model and (2) everywhere for 
the eddy viscosity model „ 


Numerical Solution of the Finite Difference Equation 

The nonlinear system of simultaneous algebraic 
equations given by Equation (3a 18) has been solved by the 
Gauss-Seidel successive substitution technique. The 
numerical scheme and the computer code developed by 
Gosman, et al . [38] , has been modified for the back-step 

problem. Successive over-relaxation or under-relaxation on 
the dependent variable, (j), in Equation (3.18) can be readily 
employed. Both relaxation techniques were used and are 
discussed later in this chapter. 
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Equation (3.18) is solved for the vorticity, 
stream function, turbulence kinetic energy, and the 
turbulence length scale in that order. Referring to Figure 
3-3, the whole flow field is swept row by row in the 
direction 5 to 4 from the upper boundary to lower boundary, 
updating values at the interior nodes using Equation 
(3.18) and at the wall nodes by the substitution formulae 
for the boundary conditions. One iteration cycle is 
complete when all the four equations have been solved and 
the <J>-values at the boundary nodes have been updated. This 
iteration process is carried on till either of the following 
two convergence criteria has been satisfied: 


i - — H ^ ~ ' 

<p max 


(3.50) 


.n .n-1 
— ) 


n-1 


0^00001 


(3.51) 


max 


max 


The second criterion is necessary because when the 
magnitude of the 4>-value at a particular node becomes very 
small, the variation, at that node can still be 

larger since it depends on the values at the surrounding 
nodes; and thus the first criterion is difficult to satisfy 
even though the rest of the field has converged. 
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Node from which the iteration 
process starts; the arrows 
show the direction in which 
the process proceeds row by row 
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Figure 3-3. Distribution of the finite-difference grid over the entire flow field 





After the convergence is achieved, the u and v 
components of the velocity are calculated. The stream 
function and the velocity profiles are plotted as the final 
results . 

Grid Distribution 

The variable grid spacing is smaller near the 
walls where large gradients in 4) are expected. Moving 
farther away from the walls, the ratio of the intervals 
between the nodes has been kept constant and equal to 1.3, 
that is, the grid spacing increases in geometric progression 
away from the walls (see Figure 3-3) . 

Selection of Constants 

The selection of the so-called constants, , 
appearing in Equation (3,9), and , and shown in 

Table 3-2, page 72, is of prime concern in the current 
turbulence research [37, 38, 44, 45, 46, 47, 48, 49] in the 
two-equation modeling of turbulence such as the TKL model or 
the (k-kil) model. In the latter model, first proposed by 
Rotta [50], the dependent variable, kil , has replaced Z in 
the TKL model because it does not diffuse at a rate propor- 
tional to . In some models, eddy viscosity, v. , is 

d Z *C 

treated as a dependent variable. These models have been 
used to predict the turbulent flow situations different from 
that of the rearward- facing step considered herein, and the 
constants have been evaluated by comparing the predictions 
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of the models with the experimental data available for 

various flow situations. For instance, Ng and Spalding [44] 

predict boundary layer flows near walls from the (k-kJt) .. 

model by using the experimental results for homogeneous 

shear flows in local equilibrium* A value for equal to 

Oel was obtained thus. The turbulent Prandtl number, , 

was taken as unity. Using the (k-kJl) model, Rodi and 

Spalding [47] obtained = 0.09 for free shear flows while 

was taken as 2*0, Launder, et al . [48], used the 

3/2 

turbulence kinetic energy dissipation rate, e = k ' /%, as 
the dependent variable instead of Jl, In the (k-e) model, 

= 0,09 and - 1.0 were found to produce satisfactory 
results for various turbulent flows [48] , To treat low 
Reynolds number flows, Jones and Launder [49] used 
= 0.09 X exp [-2 . 5/ (l+R^/50) 1 and = 1.0. 

A three-equation model [48] was used where the 
third dependent variable considered was the turbulent shear 
stress, pu' V In this model, which was also used to 
predict two-dimensional axisymmetric jets and wakes, was 
taken as 0*09. Wolfshtein [45] used C = 0.22, C. = 0.416, 

r* ^ 

and a, = 1.53 for the Couette flow case, 
k 

The incompressible, two-dimensional wake flow is 
one of the most difficult to predict properly [46], 

Apparently, the main difficulty is in predicting the 
asymptotic rate of decay of the wake, which is a weak-shear 
problem. Due to the paucity of knowledge regarding these 
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"constants" for recirculating flows, resort has been taken 
to the previous information and the computational experi- 
ments which are conducted in this study to furnish a 
reasonable set of constants for the TKL model. Unfortu- 
nately, however, these constants are not only different for 
different flow situations , but vary from one region to the 
other in the same flow situation as a function of when 
is small. As becomes very large, and 

take on constant asymptotic values, provided the hypothesis 
that k and i are adequate to characterize the state of 
turbulence is valid. However, the o's are normally, as in 
this study also, taken as unity [38]., The following 
relations define the various constants used for the present 
investigation . 


where 


where 


Cb = s ^ ' 

oo 


(3.52) 



1,0 , 



+ 1/R^ , 


(3.53) 



1.0 . 


C 


s 


Cg + 1/R^ , 

OO 


(3.54) 
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where 


C = 0,35 . 

Corresponding to Equation (3.12a), 


The above set of constants seems to be a reasonably 
good selection to begin with. Further discussion on these 
constants will appear in the parametric study reported later 
in this study. Efforts to establish some kind of a 
universal set of constants for recirculating flows are still 
being carried on. 

Discussion on Convergence, Accuracy, and Economy of 
the Solution Procedure 

Convergence . It has been experienced [38] that if 
a nonuniform grid spacing is used near the walls, both the 
vorticity b.c, (Equations (3,29) and (3.30)), when used 
explicitly for updating the vorticity at the wall, may cause 
divergence due to the coupling of the vorticity and the 
stream function equations through the vorticity b.c. 

Gosman, et al n [38], suggested as remedial measures (1) to 
hold the ratio of the intervals between the nodes (Figure 
3-3, page 98) less than 1.5 especially near the wall, where 
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this ratio should be as close to unity as possible and (2) 
to use the modified substitution formula for (o (Equation 
(3.31)). 

In calculations carried out in this investigation, 
convergence seemed to be impaired when a coarse grid with 
the ratio of the intervals between the nodes was 2 and above, 
especially in the z~direction. 

The coefficients, C*s and D, in Equation (3.18) 
vary from one iteration to the next, and hence this is a 
nonlinear algebraic equation. However, from experience, 
Gosman, et al . , point out that the convergence criteria for 
a linear set of algebraic equations to converge are often 
sufficient for equations such as Equation (3-18) also. As 
there is no such criteria for the nonlinear equations, the 
convergence of the solution procedure can be based on the 
following three criteria for a linear system of algebraic 
equations : 

I I ^ I P 1. ^ every grid node. 

2. £|c|p<l,onat least one node. 

3. The variation in the C's and D from one 
iteration to the next is small. 

The third criterion restricts the nonlinearity of 
the coefficients to an extent that they could be treated as 
linear in the whole iteration procedure, thus making the 
first two criteria meaningful. 
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It can be seen that ^ | C | p = 1 on every node in the 
field if c = 1 (Table 3-1, page 65) , which is true for all 
but the vorticity equation where c = • For low 

Reynolds number R, the value of at the node P is close 

to those at the surrounding nodes because the gradients are 
small (refer to the section on discussion of results) , and 
hence J|c|p is close to unity. But for higher R when the 
gradients in velocity are large, especially in the 
recirculating region, this criterion fails to be satisfied 
for the vorticity equation when using the eddy viscosity 
model (Equation (3.8)) because of the large variations in 
^eff neighboring nodes. That is why in the 

computation at high R, using the eddy viscosity model to 
describe turbulence, divergence was encountered. The remedy 
for this lies in under-relaxing the variable, o), or the C's 
and D coefficients in the vorticity equation. (Further 
discussion on this will appear later.) To carry on the 
computation at higher R, the TKL model was used for which no 
such difficulty was encountered. 

If the central-difference schemes for the first 
order convective terms were used instead of the upwind 
differencing, the coefficient could become negative [38]; 
and this could cause to become negative, especially at 

high R when the stream function values are correspondingly 
large and hence |a„| > cB„ . Owing to this fact, the 

xh hi 
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central-difference schemes fail to achieve convergence at 
high R due to their inability to satisfy the first 
convergence criterion,, 

Comparing a flux- type b.c., (J) = constant, with 
Equation (3ol8) , it can be seen that hence the 

second criterion is satisfied for all but the o) equation, in 
which case no flux-type boCo is used anywhere on the 
boundaries except for the case when a zero vorticity is 
used at the sharp corner 2, corresponding to the b.c. given 
by Equation (3.30). However, a comparison of the solution 
corresponding to the b.c. given by Equation (3.29) to the 
solutions corresponding to the discontinuous b.c. and zero 
b.Co given by Equation (3.30) showed no marked change in the 
rate of convergence o This indicates that either the second 
criterion is satisfied somewhere at an interior node in the 
field or else the criterion for Equation (3.18) is generally 
sufficient and not necessary for convergence; otherwise, 
divergence should have been encountered. 

The third criterion can be satisfied by under- 
relaxing the dependent variable, (j). As a better alternative, 
the C's and D coefficients should be under- relaxed as they 
are functions of quantities like as well as (j). In the 

present computations with the eddy viscosity model 
incorporated at high values of R, large variations in the 
source term for the stream function equation (Table 3-2, 
page 72) were found to occur, and hence the third criterion 
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was not satisfied. This was detected from the solutions 
which diverged. Although the strong variations in the 
source term were first seen as the cause for divergence, it 
is interesting to note that the failure to satisfy the first 
and the third criterion is really coupled through the strong 
velocity gradients and hence through strong variations in 

^eff 

The source terms in the k and H equations were 
under-relaxed in the sense of reducing the iteration-to- 
iteration variations by recasting Equation (3,18) into 
Equations (3.18a) and (3,18b), respectively. 

In the original form of Equation (3.18) for k, 
the 3/2 power of kp induces large variations in the source 
term and could cause divergence. In the present computa- 
tions, using Equation (3.18) for Z, it was observed that Z 
became negative early in the iteration process due to the 
comparatively large magnitude of the term, 


V 


P 


{ 2,W 


sktS^ 


/k Z AB , 
P 


which in turn made k negative in Equation (3,18a). To stop 
£ and k from becoming negative, Equation (3,18) for Z was 
rearranged to give Equation (3.18b) ^ Also, in Equations 
(3.18a) and (3.18b) the denominator is always greater than 
unity, and hence the first and second criteria are 
unconditionally satisfied. The third criterion is satisfied 
as the source term variations are small for the Z and k 
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equations. Obedience to these convergence criteria resulted 
in a very rapid rate of convergence of the solution for the 
Z and k equations, which testifies to the meaningfulness of 
the above convergence criteria. 

Accuracy and economy . The accuracy and economy of 
the solution procedure are opposed to one another. The 
finer the grid size and hence the more the computing time, 
the better the accuracy or the lesser the truncation error. 
Wherever one-sided difference, that is, first-order 
approximation to the first order derivatives has been used, 
steep gradients are "smeared" especially near the walls. 
Therefore, it is necessary to use a very fine grid size near 
the walls and wherever else the gradients are expected to be 
large. Finer mesh size also satisfies the assumption of 
linear vorticity near the walls. Wolfshtein [45] found by 
comparing the exact and the finite-difference solutions fv;r 
simple cases of Couette flow, impinging jet flow, and 
uniform velocity flow that a finer grid away from the walls 
does not make the vorticity solution any better; but the 
stream function in the middle of the flow field, perhaps due 
to the strong nonlinearity of its distribution, is sensitive 
to the grid size. However, by experience it was found that 
the stream function solution behaved better than the 
vorticity solution in general; and hence a reasonable 
distance away from the walls , coarser grids were used to cut 
down the computing time. 
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The effect of smearing is to introduce an addi- 
tional "false" diffusion of cj) into Equation (3.18), which 
can be assumed to be represented approximately by the 
magnitude of the "false exchange coefficient" by the 
following relation for all types of flows [45]: 


r 

r 


false 

eff 


0.36 R 


eff 


°eff ^ sin(2a) 


(3.56) 


where An is the square mesh size; a is the angle which the 

streamlines make with the coordinate system (the mesh lines) ; 

is local effective Reynolds number based 

on the step height, h; and - ^eff^^eff'^® effective 

Prandtl or Schmidt number. 

As the smearing is associated with the first order 

derivatives in the convective terms , it becomes important 

when the viscosity is very small. Hence Equation (3.56) .is 

very demonstrative of how the false diffusion effect can 

influence the solution in laminar flows as u , An, or a 

r 

(from zero to 45°) increase. For laminar flows Equation 
(3.56) reduces to 


Ffaise ^ 0o36 pu^An sin (2a) . 


(3.56a) 


Defining cell Reynolds number, R = pu An/p, Equation 

O 6 ^ IT 

(3„56a) becomes 


ffaise ' ^cell ^ sin(2a) , (3.56b) 
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Therefore, for laminar flows the importance of the two 
important parameters, ^^ell becomes clear. In the 

literature, for laminar flows the cell Reynolds number has 
been used as a parameter to control the accuracy 
(truncation error) of the method used. As a increases from 
0 to 45®, the false diffusion effect increases from zero to 
a maximum, which shows that the curvature of the streamlines 
is responsible for magnifying this effect. Thus for higher 
values of R the remedy lies in (1) choosing the correct 
streamline coordinate, a = 0°, and (2) reducing the mesh 
size, An. 


From (2) it can be seen that a given accuracy, 
that is, the truncation error within a reasonable tolerance, 
can be achieved at the cost of the computing time. 

For turbulent flows, however, as is usually 

much larger than F (for example, for the co equation, 

> u) , and is close to unity, false diffusion is 

not very important. 

The criterion for convergence has been determined 
from experience by observing that the solution obtained 
after a definite number of iterations remains essentially 
the same for further iterations. It has been experienced 
that the number of iterations required for convergence 
increases with an increase in the number of grid nodes; and 
hence, contrary to what one might expect, the computing time 
increases more than proportionally to the number of grid 
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nodes. The number of grid nodes should therefore be kept as 
low as possible 3 It is a good practice to provide the 
initial conditions discreetly in order to quicken the rate 
of convergence. It has also been experienced that the flux- 
type boundary conditions are better than the gradient-type 
boundary conditions as the latter slacken the rate of 
convergence. So wherever possible the flux-type boundary 
conditions have been prescribed. 

Discussion of Results 

The results obtained in the present study corre- 
spond to two models of turbulence, the eddy viscosity model 
(Equation (3.8)) and the TKL model (Equation (3.9)). In the 
former,, a parametric study of Reynolds number was carried 
outj and in the latter, the influence of variations in the 
constants defined by Equations (3.12a), (3.52), (3.53), and 

(3,54) was investigated. While in the eddy viscosity model 
the empiricism involved is that of the mixing length 
distribution throughout the flow field, the TKL model is 
dependent on the value of the constants chosen in the 
modeled terms. It has been noted [37] that the mixing- 
length theorem is not adequate to define the turbulence 
associated with the recirculating flows and that the two- 
equation turbulence model, like the TKL model, is better 
suited for such a flow situation. That this is true has 
been confirmed in this investigation. 
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With Reynolds number, R, as a parameter the 
results corresponding to the laminar solutions are discussed 
first. The results corresponding to the turbulent case with 
the eddy viscosity model introduced into the governing 
equations are discussed next. This is followed by a 
discussion of the results of the TKL model. 

Laminar solution (^Reynolds number as a parameter) . 
Results for laminar flows were achieved for step-height 
Reynolds number, u^h/v, equal to 6«9, 69, and 6,900 corre- 
sponding to the step height, h = 1 cm. The outlet stream 
function boundary condition given by Equation (3,22) and the 
implicit substitution formula (Equation (3.31)) for 
vorticity were used in these computations = Figures 3-4, 

3-5, and 3-6 are the streamline plots.. On comparing these 
figures it is seen that the recirculation bubble length 
becomes larger as the Reynolds number is increased, which is 
realistic From Figures .1-7, 3-8, and 3-9, which show the 
direction of flow (indicated by arrows), it can be seen that 
there are two contra-rotating eddies in the cavity zone. 

The corner eddy is smaller than the recirculation eddy. For 
R = 69, the length of the corner eddy is about 2.5 step 
heights along the ground and the recirculation eddy is 
spread over a length of seven step heights. One can observe 
from the data of Tam, et al . 18], that for a recirculation 

length of seven step heights, a region of negative surface 
pressure gradient which corresponds to the corner eddy was 
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Figure 3-4. Streamlines over a back step in laminar flow 
(Reynolds nuicdser = 6*9) • 
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Figure 3-5. Streamlines over a back step in laminar flow 
(Reynolds number = 69) . 
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Figure 3-6. Streamlines over a back step in laminar flow 
(Reynolds number = 6.9x10^). 
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Figure 3-7. Flow direction over a back step in laminar flow 
(Reynolds number = 6.9). 
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Figure 3-9. Flow direction over a back step in laminar flow 
(Reynolds number = 6.9x10^) . 


found to be about 2^.5 step heights long (see Chapter IV), 
Thus, although the results for such a low Reynolds number 
laminar flow cannot be compared with these experimental 
findings which correspond to a turbulent flow, one can 
observe that the ratio of the length scales of the two 
eddies is in agreement with that obtained from the experi- 
mental data, which gives one a qualitative feel of the 
relative size of the corner eddy with respect to that of the 
recirculation eddy. 

Also, with the decrease in R the flow seems to 
separate from a lower location down the base, which is 
meaningful because at R 0 the flow would not separate and 
remain attached to the base right down to the ground. This 
shift of the separation point down the back side of the step 
has also been found computationally by Roache and Mueller 
[41] „ As it is not known from what point the flow separates 
on the base, the plot of the zero streamline is started from 
the interior point immediately next to the base. At the 
inlet for R = 6.9. and 69, the upshoot in the streamlines can 
be attributed to the fact that the logarithmic velocity 
profile used is not a realistic boundary condition for such 
low Reynolds nximbers ^ But for R = 6,900, the inflow b.c. 
is more meaningful as can be seen from Figure 3-6. 

It may be pointed out that at R = 6,900 the 
central-difference schemes used for the first order 
derivatives in the x-direction will render the method 
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unstable as discussed earlier. Hung and Macagno [52] , in 
spite of the use of unsteady equations which are believed to 
be more stable than the steady ones, could not obtain 
results for. Reynolds numbers higher than 333 (using the 
three-point central-differencing formula in the analysis). 
With the upwind-differencing used in this study laminar 
solutions for R still higher than 6,900 could have been 
achieved. The present study, being primarily concerned with 
turbulent flow, was, however, at this point directed towards 
turbulent solutions. Toward this goal the eddy viscosity 
model was introduced into the governing equations . 

Eddy viscosity model (Reynolds number as a 
parameter) . For Reynolds number equal to 6.9, the eddy 
viscosity model yielded a "turbulent" solution which was 
more or less like the laminar counterpart. The reason for 
this can be explained as follows: The maximum velocity 

gradient, , at the ground, z = 0, for a surface 
roughness, z^ = 0^01 cm corresponding to the step height, 
h = 1 cm, is of the order of 2.0 (see Equation (3.57)). 
Intuitively it can be seen that the order of magnitude of 
(■^) at the base should not be very much different from 

C7 

3 Ul 

that of (-^) at the ground. Thus it follows from 

d z max 

Equation (3.8) that << pi, which explains why the solution 
remains practically unchanged with the introduction of the 
eddy viscosity model. 
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For the turbulent flow calculations the laminar 
solution was taken as the initial conditions which resulted 
in rapid convergence « Turbulent solutions obtained, 
beginning with the same initial conditions as used for the 
laminar solutions, took much longer to converge. 

Next, the eddy viscosity model was introduced into 
the laminar solution for R = 6,900. The turbulent solution 
in this case diverged. The cause for divergence was traced 
to the following; 

Considering u^ as the reference velocity at the 
upper right-hand corner of the control volume (point 5 in 
Figure 3-2, page 82), Equations (3.4) and (3.38) give 


3u _ 


u 


(z+z, 


r 


)£n 





(3.57) 


which shows that the velocity gradient, , will be steeper 
for higher values of u^ ; thus from Equation (3.8) will 
be large and correspondingly its variation from node to node 
will be large. Thus, as discussed earlier in the section on 
convergence of the solution procedure, divergence was 
encountered in the vorticity equation and eventually in the 
stream function equation also. 

To overcome this difficulty two remedies are 

possible : 
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1 


o Reducing the mesh size such that the u 

erf 

variation from node to node is sufficiently 
small. 

2. Using under- relaxation on the C's and D 

coefficients in the vorticity equation or on 
the vorticity, w, itself. 

As the criterion for selecting a proper mesh size 
is not known, it is easier to make use of the under- 
relaxation technique. Under-relaxation was tried on the 
03-equation for the case of a higher Reynolds number, 

g 

R = 6,9x10 , corresponding to the step height, h =1 meter. 
The under-relaxation parameter, p = 0.01, was found suitable 
to use after experimenting with different values r With 
p = 0,01, vorticity, o3, was calculated at the end of each 
iteration as follows: 

n . n-1 . 

03 — p03 + 03 (1-p) . 

The turbulent solution for R - 6,9x10^ was thus obtained 
(Figure 3-10) , In this solution it was observed that 
although the rest of the flow field had converged according 
to the convergence criterion of Equation (3.50), at the 
point half-way down and next to the base the value of 
vorticity was vanishingly small, and hence the above 
convergence criterion was not acceptable. Instead, the 
convergence criterion given by Equation (3,51) was used as 
discussed earlier in this chapter. The solution gave a 
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Logarithmic velocity profile, 



at the outlet 



Figure 3-10. Streamlines over a back step in turbulent flow using 
the eddy viscosity model (Reynolds number = 6.9x10^). 


recirculation length which is too large.. This is due to the 
eddy viscosity model being inadequate to correctly predict 
the turbulent recirculating flows « Thus a better model was 
sought. 

Over- relaxation on the stream function was also 
used [38] , but it did not seem to produce appreciable 
increase in convergence rate, 

TKL model . In this model the values of the 
constants, C » discussed earlier, 

must be determined by comparing the predicted results of the 
present model with the experimentally known results. 
Experimental results such as the length of the recirculation 
zone, the location of the dividing streamline, and the 
variation of turbulence intensity in the cavity zone, 
obtained for a back step in a wind tunnel for laboratory 
flows and for wind towers in atmospheric flows are the 
measured quantities of interest for comparison purposes. 

Referring to the modified substitution formulae 
for k and it, Equations (3.18a) and (3.18b), an increase in 
in the "dissipation" term reduces the turbulence kinetic 

^00 

energy, k, which in turn induces a decrease in the length 

scale, £, Thus from Equation (3,9), it is observed that 

VI decreases with an increase in , Similarly, an 

00 

increase in C, appearing in the "eddy breaking" term and a 

00 

decrease in C in the "eddy stretching" term of the length- 

®oo 

scale equation induce a decrease in , which means the 
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turbulence is decreased 9 This can be interpreted as the 
tendency of the flow to behave more like a laminar flowo It 
has already been observed that the higher the R, the larger 
the bubble in laminar flow. Thus changes in C, , C , and 

00 ®oo 

C, can be made to control the size of the recirculation 
bubble y and also the level of turbulence kinetic energy. 

The effect of C , on the other hand, is not so obvious. 

^00 

It would seem on a casual inspection of Equation (3.9) that 

an increase ii' C will increase proportionally, but 

this is not the case. In the present problem the turbulent 

Reynolds number, R. >> 1, in general, hence C - c . From 

^ ^00 

Equations (3c,43), (3.44), (3.48), and (3.49) it can be seen 
that on all the boundaries an increase in C results in a 
decrease in the value of k by the factor, l/C^ ^ Looking at 

^OO 

Equation (3ol8a), it is observed that k values at the 

tr 

interior nodes will not decrease as much as the k values at 

1/2 

the boundaries because of the k ' term in the denominator.. 
From Equation (3.18b) Z also decreases slowly compared with 
the increase in C . Therefore, from Equation (3.9), 

^oo 

although the value of on the boundaries does not change 

with C , the value of at the interior nodes will 

^00 

increase slowly; this has been found in the computation 
carried out in the present study. 

Knowing approximately the influence of these 
constants on the solution, the next step is to vary them 
judiciously to achieve results in agreement with experi- 
mental data. 
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The boundary conditions on k at the outlet. 

Equations (3,44) and (3.44a), have both been investigated. 

The latter folating b.c. with the corresponding boundary 

condition for ip at the outlet. Equation (3.23), did not seem 

to affect the rate of convergence over that of the former, 

and it yielded more parallel and horizontal streeunlines ; 

therefore, it was used for most of the results shown. The 

reason the gradient-type b.c. did not slacken the rate of 

convergence is that the outlet has been taken sufficiently 

far downstream so that the gradient effects are negligible. 

At the upper boundary a floating (gradient-type) b.c. on 
3 Z 

Z, = 0, was tried, which worsened the rate of convergence, 

and hence was dropped in favor of the fixed b.c. there. 

Initially, all the constants were set equal to 

unity. This gave an unstable solution with the stream 

function becoming negative over most of the flow field. 

Different values of the constants were tried, based on the 

discussion of their influence on the k and Z equations. 

Figures 3-11, 3-12, 3-13, 3-14, 3-15, and 3-16 show the 

stream function plots for various combinations of these 

constants. It can be seen that for the range of values of 

C chosen (0.05 to 1.0), the recirculation bubble does not 
^00 

change appreciably. On the other hand, the constant, C , 

®00 

appearing in the "stretching" term for the length scale 
equation, has a very marked influence on the size of the 
bubble. For =0.5 the recirculation zone length is 
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Logarithmic velocity profile, u 
specified at the outlet 


t 





in turbulent flow 
6.9x106) 
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2+Zq 

Logarithmic velocity profile, u = — i^n t 

specified at the outlet 
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Logarithmic velocity profile, u 
specified at the outlet 


K “ Zq 


4 


Floating boundary conditions used at the outlet: 
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Figure 3-14. Streamlines over a back step in turbulent flow 
using the TKL model (Reynolds number = 6.9x10®), 

C = 0.05, C, = C , = 0.1, C = 1.0, 

^00 00 00 oo 


Floating boundary conditions used at the outlet: 



Figure 3-15, Streamlines over a back step in turbulent flow 
using the TKL model (Reynolds number * 6.9x10 ), 

"oo *^00 OO OO 
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condl'tlons used at the outlet: 
0 



Floating boundary 


3^i|) 


rt • J.U 

= 0 WXth -gj 



0 • 4 the 


about 2el step heights (Figure 3-15), for C = 

®oo 

length is 2.5 step heights (Figure 3-16), and for C = 0.35 

®oo 

it increases to 4.2 step heights (Figure 3-17). In the 
computations carried out with lower values of C than 0.35 

®oo 

it was observed that the boundary condition on k at the 

outlet (Equation (3.44a)) tended to impair the stability of 

3k 

the solution; and hence an alternative b.c. , = 0, was 

used at the outlet to investigate the effect of further 
decrease in C on the solution, which appreciably increased 

®oo 

the size of the recirculation zone downstream. 

The value of 0,1 for C was found more suitable 

^ 00 

than the values of 0^08 and 0.05 as far as the rate of 
convergence of k and i equations is concerned. Also, this 

value of C has been used by other authors in various 

u 

turbulent flow situations (see the section on selection of 
constants) o The value of loO for and Cj^ has been 

00 <x> 

chosen out of the experience gained by carrying on the 
computational experiments. Thus with the following set of 
constants , 


and 


C - 0.1 , 


C = 0.35 , 
s_ 


^b C - 1.0 , 

OO 00 
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Floating boundary conditions used at the outlet: 



Figure 3-17. 
using the 
C = 0.1 
^00 


I 


streamlines over a back step in turbulent flow 
TKL model (Reynolds number = 6.9x10”), 




= 1 . 0 , 



0.35. 


the recirculating flow behind a back step has been well 
predicted; and they could be used for recirculating flows in 
different flow situations also. These asymptotic values for 
the functions, ' ^d ' and , as may be the 

"universal" constants for recirculating flows. 

Conclusions and Summary 

The present approach of solving the two-dimensional, 
incompressible turbulent Navier-Stokes equations in their 
steady form as applicable to atmospheric boundary layer flow 
has been shown to yield results for recirculating flows 
behind a backward- facing step in agreement with the experi- 
mental findings. Whereas for low and intermediate Reynolds 
numbers the eddy viscosity model of turbulence introduced 
into the governing equations has proved satisfactory, it has 
failed to predict the real flow situation at high Reynolds 
numbers when the magnitude of the velocity gradients in the 
flow field is correspondingly large. This is attributed to 

the fact that in recirculating flows the relation between 
stresses and the velocity gradients is unknown and more 
complicated than that assumed by the eddy viscosity model, 
and that where the velocity gradients vanish in recirculating 
flows, turbulence does not necessarily become zero* On the 
other hand, the TKL model of turbulence has been demon- 
strated to produce realistic results at high Reynolds 
numbers. The influence of the coefficients in the modeled 
source terms in the turbulence kinetic energy and turbulence 
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length scale equations on the predicted results has been 
computationally investigated^ and a set of constants has 
thus been found for the problem under study o 

It is felt that even more accurate predictions can 
be made with a three-equation model of turbulence in which 
the third dependent variable can be taken as the turbulent 
shear stress, pu' v^ . 

For laminar flows the recirculation bubble was 
seen to increase in length with Reynolds number, which is 
physically trueo The two contra-rotating eddies resulted in 
the solution, the corner eddy being the smaller one and 
rotating in the opposite sense to that of the larger one, 
the recirculation eddyo 


CHAPTER IV 


APPROXIMATE INTEGRAL TECHNIQUE 
Empirical Information 

The approximate integral method to predict the 
velocity profile in the recirculation region behind a 
rearward- facing step is based on the experimental infor- 
mation given in [4, 8, 10, 21J , The available data 
correspond to a "strong perturbation" (see Chapter I) and 
have been obtained from experiments conducted in wind 
tunnels at controlled turbulence levels. The results from 
these investigations correlate well, and the present 
mathematical model has been shown to reproduce the different 
features of the flow. In addition, a surface eddy viscosity 
distribution has resulted in the model which follows the 
empirically known surface pressure gradient distribution 
curve.. Figure 4-1 shows the surface pressure coefficient 
plotted against the nondimensional distance, x/h [8], where 
h is the height of the step, is the velocity in the 
undisturbed flow, 6 q is the boundary layer thickness at the 
separation point, and x^^ is the reattachment length. 

Figure 2-6, page 13, shows the location of the reattachment 
point (CR = x„) and the location of the dividing streamline * 

The error function profile which is known to 
approximate the fully developed shear layer profile in the 
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Figure 4-1. Varia 
coefficient be] 


a^i 


mixing region is given by 


u 


u 


erf ~ ”2” 


oz. 


(1 + erf 


where is the undisturbed velocity which is asymptotically 

o z 

approached as the transverse coordinate, — and a is 

the spreading rate parameter; also, 0 as ^ 


Governing Boundary Conditions for the Problem 

To develop the governing equations consider a two- 
dimensional, steady, incompressible turbulent flow over a 
rearward- facing step, ABC (Figure 4-2) . The separated flow 
is modeled as shown. 

The initial boundary layer at the separation point, 
B, is considered to be represented by the simple power law 
profile , 




n 


where u^ is a reference velocity, the index n = 7, and 6^ is 
the initial boundary layer thickness. As the flow separates 
at B, a new shear layer originating there develops and 
spreads linearly downstream so that its upper boundary is 
represented by BB' - The dividing streamline BER coincides 
with Bx' up to some point, E (refer to Figures 2-3, page 10, 
2-6, page 13, and 2-8, page 23). BER encloses the sepa- 
ration bubble in which there are two contra-rotating eddies , 
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the small corner eddy and the large recirculation eddy. DD* 
represents the boundary of the slowly growing initial 
boundary layer of almost constant height (see Figures 2-3, 
page 10, and 2-21, page 48) . Camarata [21] has successfully 
correlated the error function profiles with the experimental 
turbulent free shear layer profiles (see Figure 2-24, page 
52) by making use of the conservation of momentum principle. 
Also, the error function profile is shown to match with the 
outer boundary layer power law profile smoothly except for 
the upstream segment of the pre-asymptotic region where the 
slope-matching is discontinuous. 

As a first approximation, Bx* is considered as the 
locus of the u/u^ = 0o5 points » Taking Cx and Cz as the 
abcissa and the ordinate of the Cartesian coordinate system, 
respectively, the velocity in the new shear layer will be 
represented by the error function profile o 


u 


erf 


u 

e 

T 


(1 + erf 


g (z-h) 

X 


} 


Referring to Figure 4-3 (EF is a control surface 
with a unit depth perpendicular to the plane of paper) , 
momentum is conserved between station x = 0 and any other 
x-station. Assuming that static pressure variations have 
negligible influence on the momentum balance and that the 
net momentum in the recirculating region is small, the 
momentum equation becomes 
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Figure 4-3. Illustration of the conservation of 
momentum principle. 


^ 2 

/ u dz 
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x=0 0 


in "Q rr / V> ^ ^ ^9 

S t-^ (1 + erf 2i|lEL}] dz + / u'^dz 
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® / {1 + erf dz . (4.1) 
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For a continuous velocity profile, 


" = "erf . at z = Zj„ , 


or 


or 


z u a(z -h) 

Ur(^) ^ ^ a ^ erf , 


1/7 


u 


2u (/i) 


1 + erf 


TT^^ET • 

m 


(4.2) 


For a given c, Equations (4,1) and (4.2) yield z^ 

and u . Thus the error function profile is completely 
e 

specif iedo 

Referring to Figure 4—4 , it is assumed that the 
recirculation zone velocity profile is expressed by a fifth 
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degree polynomial as follows : 

u _ = e(x)z + f (x)z^ + g(x)z^ + h(x)z^ + k(x)z^ , (4.3) 

3TZ 


where for a given x the coefficients, e, f, g, h, and k, are 
constant. This profile should match with the new shear 
layer profile at a point, z^ , on the dividing streamline. 
From experience it was found that for a continuous matching, 
the following three conditions should be satisfied at z = 

^rz '^erf ' 


2 . 


3. 


3u 


rz 

Tz~ 


au 


erf 

IT' 
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a^u 


rz 
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u 


az 


erf 

• 


Mechanics of the flow in the recirculation regidn 
suggests that the equation of continuity should be satisfied 
there which provides the following condition; 


4. / ^ u dz = 0 0 

rz 

0 

Also, Equation (4.3) automatically satisfies the no-slip 
condition at the ground, z = 0® 

Introducing the x-direction Navier- Stokes equation 
(Equation (3.2)) at the ground gives 
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i ^ sr -i. 
p dx 3z 


Su. 


[u 


rz 


eff 


where Pgff is given by Equations (3n5), (3.6), and (3.7), 

From Equation (4.3)atz=0, 


3u 

rz 

"^z~‘ 


= e 


Therefore, at a given x- location. Equations (3.6) and (3.7) 
yield the following expression for eddy viscosity; 


V 


t 



(4.4) 


and thus the momentum equation yields an additional 
condition in the following form: 


5. 


i ^ 

p dx 


- {kz, 


3 u 


+ v) 


3z 


rz 


The aforementioned five conditions along with 
Equations (4c 1) and (4,2) constitute the governing boundary 
conditions for the solution of the recirculation zone 
velocity profile, and upon simplification take the following 
form, respectively. 

eZj^ + fZj^^ : + 9 ^ 1 ^ + ^^ 2 .^ " '^erfl ' (4.5) 

|Zl 
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where 


u 


erf 


u a(z,-h) 

^ u + erf 


} , 


4 ®“erf 


e + 2fz, + 3gz,^ + 4hz,^ + 5kZj^’ = — S| 


/ (4.6) 


where 


3u 


erf 


Tz 


x/tt 


o{z.-h) 

[exp - { — ] , 


:)2 

3 u 


2f + 6gz^ + 12hz^ + 20kZj^ 


_er f 
dz 


(4.7) 


where 


9^u 


dz 


erf 

T 


2a (z.,-h)u 
1 e 


o (z,-h) 

[exp - { } ] , 


30e + 20fz^ + 15gz^^ + 12hz^^ + lOkz^^^ = 0 , (4.8) 


f = 



(4.9) 


where 


4 {v^+v} 



Solution of Equations (4.5) through (4.8) leads to the 
expression for the coefficients of Equation (4.3) given in 
Table 4-l« 


Normalized Foimi of the Equations to be Solved 

Nondimens ionali zing the velocities « u, u « u. f 

JTZ 6 

and Ugyf i with respect to the reference velocity# u^, , and 
the distances# x# x^^ # z# z^ # and # with respect to the 
step height# h# gives 
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Therefore, Equations (4,1), (4^2), and (4.3), respectively, 
take the following forms: 


0(2,-l) ^ 

— } dz* , (4,1a) 
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TABLE 4-1 


EXPRESSIONS FOR THE COEFFICIENTS OF EQUATION (4.3) 




/<Sn ) 


1/7 


u 


1 + erf 


X* 


(4.2a) 


and 


u 


rz, 


e*(x*)z* + f*(x*)z* + g*(x*)z* + h* (x.) z* 


+ k* (X*) 2* . (4.3a) 


Equation (4,4) in the nondimens ion a 1 form is 


e*| = 


(4.4a) 


Equation (4.9) becomes 


hUj.dCp/dx* 

f* 


(4.9a) 


The expressions for the coefficients, e* f g* # h* 
and k* , remain the same as given in Table 4-1 except that 
all the quantities used therein are nondimens ional , 

Thus the problem is reduced to solving Equations 
(4.1a), (4p2a), and (4p3a) with Equations (4.4a) and (4.9a) 
and the coefficients of Table 4-1 in the nondimensional form 


Correlation of the Spreading Rate Parameter with the 
Experimental Data 

The error function velocity profile which would 
approximate the free shear layer velocity profile is assumed 


149 



to tend to minus infinity in the negative z-direction [ 2 , 3 , 
21 , 23]. This imposes a restrictive condition on the 
selection of a as the transverse coordinate in the error 
function velocity profile should approximately follow the 
relation. 


at the ground, z* 
that is, 

Also, as a determines the rate of spreading of the new shear 
layer, the development of the line BB' in Figure 4-3, page 
141, depends upon the value of o chosen Figure 2-3, page 
10, shows that the upper boundary of the new shear layer has 
a dip near the region of sudden pressure rise (Figure 4-1, 
page 137), that is, the region of rapid distortion as 
defined by Bradshaw and Wong [4] . The line BB* thus has a 
downward- going trend until the reattachment point, when the 
effects of the distortion die out; and it rises again into 
the initial boundary layer downstream. 

It may be noted that a smaller a gives a larger 
spreading, and vice versa. Thus either a constant value or 
a very slowly increasing or decreasing value of o may be 
used up to the region of maximum pressure rise subject to 
the condition of Equation (4ol0). Beyond this point, an 


> -2 , 


= 0, and for a given [erf(-2)~ -1.0], 


0 ^ 2x^ . 


(4.10) 
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increasing value of a subject to the condition of Equation 
(4.10) should be used. In Chapter II it has been discussed 
that varying values of a have been used for the case of 
back steps in incompressible flow (range of variation from 
9 to 15 for a) • 

Solution Procedure 

With this information available for the spreading 

rate parameter, a, one can proceed to solve the system of 

the integral equation (Equation (4.1a)) and Equations (4.2a), 

(4.3a), (4.4a), (4.9a), and those given in Table 4-1, page 

148, in the nondimensional form. 

For a given a, Equation (4.1a) with Equation 

(4.2a) is solved using Simpson's fifth order quadrature 

formula for integration. The value of is calculated 

m* 

iteratively. The solution of these two equations gives 

z and u , and hence the error function velocity profile 
m^ 

is completely known. The coefficients in Table 4-1 in the 
nondimensional form are calculated next and then the 
recirculation zone velocity profile is determined by 
Equation (4.3a). 

The coefficient, e* , as calculated from Table 4-1 
is compared with Equation (4.4a) as follows: Denoting e* 

from Table 4-1 by El, for compatibility the following should 
be true : 
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|El| = 0 , (4,11) 

(KZq^ h u^) 

The nonlinear equation (Equation (4 oil)) is a monotonically 
decreasing function of , and it can be solved for 
iteratively o A simple secant method is sufficient to find 
the zero of this function. Thus the solution results in the 
value of the eddy viscosity o 

The computer code which has been developed to 
solve the present problem is explained by means of a flow 
chart as shown in the appendix. 

Comments on Convergence, Accuracy and Economy of the 
Integral Method 

The approximate technique as discussed above is a 
very fast method of predicting the velocity profiles and the 
surface eddy viscosity distribution in the recirculation 
zone behind a back step. The average time for a single 
x-station computation on an IBM 360/65 computer is about 
two seconds o 

Equation (4.1a) has two solutions for z , as 

m* 

shown in Figure 2-22, page 49. The matching is discontinuous 

at the lower value of z j but at the higher value, which 

m* 

is taken as the true solution, it is smooth. However, in 
the region close to the base of the step, the matching is 
not smooth (see Figure 4-5) . This has been experimentally 
verified [21] , and the reason could be ascribed to the 
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Figure 4-5. i: 
error func 
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interference from the base. Farther downstream, the 
matching is smooth. Due to the non-uniqueness of the 
solution for Equation (4.1a) no standard iteration technique, 
for example, Newton's method, was used; and hence the 
solutions were determined by regularly incrementing the z* 
value from 1.0 upwards. The incremental value chosen was 
such that the magnitude of the function given by Equation 
(4.1a) at the root, z , was less than 0.001. This is a 
satisfactory convergence criterion as is evident from the 
smooth matching between the error function and the power law 
profiles shown in Figure 4-6. The convergence of the 
solution of Equation (4.11) was based on the relative error 
criterion. 


- V 


< 0.001 


where and are the two successive eddy viscosity 
values at any stage in the iteration process. 


Discussion of Results 

Referring to Figure 4-7 (a) and (b) , various values 
of a were used in the computation discussed in the section 
on its correlation with the experimental data. An eddy 
viscosity distribution along the surface was faired as shown 
in (b) , which follows the magnitude of the surface pressure 
gradient distribution shown in (a) . A constant value of o 
equal to 8 was used up to x^ = 4.0, and a «= 2x^ was used up 
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(a) Magnitude of surface pressure gradient distribution 




to X* = 6.9. This variation of a gives the spreading in the 
region of rapid distortion as shown in Figure 2-3, page 10. 
Also, the eddy viscosity, , follows 



very realistically. The peak values of in the region of 



defined as the reattachment zone, confirms the high 
turbulence and hence the high heat transfer rates in this 
zone. Near the base and very near the reattachment point, 
however, does not follow the magnitude of the pressure 
gradient, which are the regions where no measurements of 
velocity or of turbulence intensity have been made (see 
Figure 2-6, page 13), and in which regions the measurements 
of shear stress are scanty and unreliable. However, it is 
well documented that the surface shear stress and hence the 
surface eddy viscosity should increase after the reattach- 
ment point. If the reattachment point is located upstream 
of the point where 
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then the eddy viscosity distribution as shown in Figure 4-7 
is quite realistic* However, due to the lack of experi- 
mental evidence, no improvements in the present model have 
been suggested as far as the eddy viscosity distribution is 
concerned. 

The velocity profiles are shown plotted at various 
x-stations in Figure 4-6, page 155, The recirculation zone 
velocity profile matches very smoothly with the error 
function velocity profile at the dividing streamline. In the 
region of the negative surface pressure gradient, the 
presence of the corner eddy with a reversed recirculating 
flow is confirmed. The large magnitude of the velocity in 
this region alludes to the presence of a very strong corner 
eddy near the basen The magnitude of the reversed velocity 
in the region of the recirculation eddy is of the order of 
0*1 u^ (Figure 4-6), and it does not exceed 0.2 u^ at any 
station. This, is in agreement with the experimental results 
of Bradshaw and Wong [4], 

At the dividing streamline, the turbulent shear 
stress, which is also the maximum shear stress in the shear 
layer, is given by 


T 

t 


2 


U 


(4.12) 
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2 

Nondimensionalizing with respect to pu^ , 


2 


(4.12a) 


The normalized shear stress, x. , has been 
plotted against x* in Figure 4-8. Although the validity of 
the momentum transfer theory and hence the Prandtl mixing- 
length theorem (Equation (4ol2)) is not quite valid in 
nonequilibrium layers, especially in the recirculation zone, 
on comparing Figure 2-7, page 16, and Figure 2-12, page 27, 
with Figure 4-8, it can be seen that the general trend in 
the value of shear stress along the dividing streamline to 
decrease downstream towards the reattachment point, as 
exhibited by Figure 4-8, is physically correct. 


Su ggestions to Improve the Model 

The fact that the static pressure behind a bluff 
plate (see Figure 2-13, page 28) varies in the transverse 
direction entails an improvement in the application of the 
z-direction momentum theorem in the case of the back step 
behind which the static pressure variations at different 
x-stations follow closely those behind a bluff plate. Thus 
Equation (4.1) will be modified as follows: 
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(4.1b) 


where pj^ and are the base and reference (undisturbed) 
pressures, respectively* 

Solution of the problem under study is strongly 
dependent upon the value of a chosen; and therefore, the 
development of line BB* in Figure 4-3, page 141 (spreading 
of the upper shear layer) , to approximate that shown in 
Figure 2-3, page 10, would require a different distribution 
of a values when Equation (4.1b) is used in place of 
Equation (4d), hence giving a more realistic surface eddy 
viscosity distribution. Equation (4.1b) also takes into 
account the interfering effect of the base for x-stations 
close to it, and thus velocity profiles in the region of 
the corner eddy could be computed to a greater reliability. 
Also, near the reattachment point where the transverse 
pressure gradient close to the surface is appreciable. 
Equation (4olb) would yield more meaningful results. 
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APPENDIX 


The numbers appearing on the left side of the 
subroutine blocks in the flow chart will be used to describe 
the function of each subroutine in sequence 

Sv±>routine CFIT uses the IBM standard sub- 
routines DATSG and DALI. DATSG rearranges the input matrix 

of X* and C values (see Figure 4-1, page 137) in an 
P 

ascending order of x*. DALI interpolates function values of 
Cp for given argument values of x* using Lagrangian interpo- 
lation with Aitken*s Scheme [54] o CFIT returns the 
interpolated vector, , to the calling program. 

Subroutine CALCUE calculates the value of u^ 
by solving Equations (4.1a) and (4o2a), It calls the sub- 
routine INTEG which uses Simpson fifth-order quadrature 
formula for integration to calculate the integral on the 
right-hand side of Equation (4«la)o The function sub- 
programs, DEF and F, are employed by INTEG and CALCUE to 
calculate u - and u* , respectively o The value of u is 
returned to the calling program. 

Subroutine CALVAR calculates the values of 


“erf* ' — ^ 


at z 


using the equations associated with Equations (4.5), (4.6), 
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and (4,7), respectively, in the nondimens ional form, and 
returns them to the calling program. 

Subroutine FUN calls the subroutine CALCOF to 
calculate the coefficients, k* , h* , g* , and e* , in that 
sequence each time FUN is called from the calling program. 
Also, it calculates the function value given by Equation 
(4.11), and returns it to the calling program. 
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Read two initial values 
of to bracket the 
root of Equation (4«ll)yl 
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No 


Call 

FUN 
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Find the root to 
desired accuracy 
using 

Secant Method 

Calculate the coef- 
ficients , e*, f*, g*, 
h*, and k«; calculate 
u^„ profile and x. 

along the dividing 
strecunline from 
Equation (4.12a) 


Trint u* , '^erf* 
and u__ profiles; 
print V. and x. 
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